plus sc10x

i.P06 and P30.PBS
check DAM signature

  1. P30.PBS and P30.LPS
    check TDT difference/shift
load dependancies

cleaning up

processed obj

GEX.seur <- readRDS("../Spp1tdt.GEX.seur.sort1006.rds")
GEX.seur
## An object of class Seurat 
## 17900 features across 23829 samples within 1 assay 
## Active assay: RNA (17900 features, 1246 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
color.FB <- ggsci::pal_igv("default")(49)[c(8,32,40,
                                            34,26,33,28,
                                            2,43,18)]
color.FBraw <-  c("#FF6C91","lightgrey",color.FB)
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB,
        pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB,
        pt.size = 0, group.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc

filtering

GEX.seur <- subset(GEX.seur, subset= preAnno %in% c("Microglia") & DoubletFinder0.05 == "Singlet")
GEX.seur
## An object of class Seurat 
## 17900 features across 22261 samples within 1 assay 
## Active assay: RNA (17900 features, 1246 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB,
        pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB,
        pt.size = 0, group.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc

subset2

GEX.seur <- subset(GEX.seur, subset= age %in% c("P30.PBS","P30.LPS"))
GEX.seur
## An object of class Seurat 
## 17900 features across 15696 samples within 1 assay 
## Active assay: RNA (17900 features, 1246 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
color.FB2 <- color.FB[c(1:3,4:7)]
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB2,
        pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB2,
        pt.size = 0, group.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "sort_clusters")

advanced

GEX.seur <- subset(GEX.seur, subset= percent.mt < 7.5 & nFeature_RNA > 1000 & nFeature_RNA < 3200 & nCount_RNA < 9000)
GEX.seur
## An object of class Seurat 
## 17900 features across 15536 samples within 1 assay 
## Active assay: RNA (17900 features, 1246 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "preAnno")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB2,
        pt.size = 0.1, group.by = "FB.info")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        #ncol = 3, 
        ncol = 2,
        cols = color.FB2,
        pt.size = 0, group.by = "FB.info")

plota <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "age") 
plotb <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "age") 
plotc <- FeatureScatter(GEX.seur, feature1 = "nCount_RNA", feature2 = "percent.rb", group.by = "age")
plota + plotb + plotc

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0.05, group.by = "sort_clusters")

VlnPlot(GEX.seur, features = c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), 
        ncol = 2,
        pt.size = 0, group.by = "sort_clusters")

par(mar=c(6,4,2,3))
sl_stat <- table(GEX.seur$FB.info)
barplot(sl_stat,ylim = c(0,3800),
        col = c("#FF6C91","lightgrey",color.FB),
        main = "Feature Barcode statistics",cex.names = 0.75, xaxt = "n")
axis(1,1:12*1.2-0.48 ,levels(GEX.seur$FB.info), las=3, cex.axis=0.85)
text(x=1:12*1.2-0.45,y=sl_stat+245,paste0(sl_stat,"\n",100*round(as.numeric(sl_stat/sum(sl_stat)),4),"%"),cex = 0.75)

subset2

re-clustering

GEX.seur <- FindVariableFeatures(GEX.seur, selection.method = "vst", nfeatures = 1800)

# Identify the 10 most highly variable genes
#top10 <- head(VariableFeatures(GEX.seur), 10)
top20 <- head(VariableFeatures(GEX.seur), 20)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(GEX.seur)
plot2 <- LabelPoints(plot = plot1, points = top20, repel = T)
plot1 + plot2

head(VariableFeatures(GEX.seur), 300)
##   [1] "Saa3"          "Cxcl10"        "Ccl4"          "Ccl3"         
##   [5] "Ccl5"          "Marco"         "Cxcl13"        "Apoe"         
##   [9] "Ch25h"         "Lcn2"          "Ccl2"          "Ccl7"         
##  [13] "Ccl12"         "Spp1"          "Rsad2"         "Xist"         
##  [17] "Il1b"          "Cd69"          "Cxcl9"         "Egr1"         
##  [21] "Acod1"         "Slc13a3"       "Cybb"          "Iigp1"        
##  [25] "Il12b"         "Cd74"          "Ifit3"         "Hist1h1b"     
##  [29] "Lpl"           "Ifit2"         "Lyz2"          "Tnf"          
##  [33] "Cp"            "Top2a"         "Egr3"          "Cd72"         
##  [37] "Ifit1"         "Fth1"          "Fpr1"          "Ly6k"         
##  [41] "Aldh1a3"       "Hba-a1"        "Gpr84"         "Ifitm3"       
##  [45] "Hba-a2"        "Wfdc21"        "Olfr889"       "Mki67"        
##  [49] "H2-Ab1"        "Hbb-bs"        "Gm50255"       "Lgals3"       
##  [53] "H2-Aa"         "Wfdc17"        "Gm26885"       "Ccr7"         
##  [57] "AW112010"      "Ifi27l2a"      "Cxcl2"         "Fn1"          
##  [61] "Ms4a6c"        "Atf3"          "Loxl2"         "Ube2c"        
##  [65] "Tnfaip2"       "Prc1"          "Dennd5b"       "Wdcp"         
##  [69] "H2-Eb1"        "Cd52"          "Gm26917"       "Trap1a"       
##  [73] "Ifi211"        "Pbld1"         "Kcnk4"         "Gdf15"        
##  [77] "Ifit3b"        "Olfr1369-ps1"  "Oasl1"         "Isg15"        
##  [81] "Gadd45b"       "Hmgb2"         "Ptges"         "Cxcl16"       
##  [85] "Egr2"          "Ccr1"          "Hmmr"          "Hsp90ab1"     
##  [89] "Cenpf"         "Ptgds"         "Fabp5"         "Ifi203"       
##  [93] "Hist1h2ap"     "Apoc1"         "Cd79a"         "Ms4a4c"       
##  [97] "Tlr2"          "Dqx1"          "Ifi204"        "Nusap1"       
## [101] "Il10"          "Hbb-bt"        "Esco2"         "Gm29773"      
## [105] "Saa1"          "E2f7"          "Slfn5"         "Cpd"          
## [109] "Kif11"         "Cdca3"         "Ifi44"         "Nr4a1"        
## [113] "Lockd"         "Bst2"          "Hspa5"         "Vpreb3"       
## [117] "Gm10457"       "Lilrb4a"       "Kcnq1ot1"      "Pclaf"        
## [121] "Gm19409"       "Ifi207"        "Csf1"          "Ifitm7"       
## [125] "Cd63"          "Tsix"          "Pbk"           "Ier3"         
## [129] "Cd83"          "Mir155hg"      "Gm8251"        "Sod2"         
## [133] "Clec4a1"       "Lhfpl4"        "Ptprcap"       "Rhox5"        
## [137] "Cdkn1a"        "Msr1"          "Dnase1l3"      "B930036N10Rik"
## [141] "Ifnb1"         "Gm34455"       "Plk1"          "C4b"          
## [145] "Gm15987"       "Stmn1"         "Srgn"          "Spib"         
## [149] "Ifi205"        "Ifi47"         "Ly6a"          "Rrm2"         
## [153] "H2-K1"         "Hist1h1a"      "Vgf"           "Slc1a2"       
## [157] "Hsp90aa1"      "Dennd3"        "A330032B11Rik" "Vnn3"         
## [161] "Oasl2"         "1700120C14Rik" "Hist1h2ae"     "Plac8"        
## [165] "1500002C15Rik" "Usp18"         "Mmp12"         "Ptn"          
## [169] "Fcgr4"         "Gm49756"       "Ms4a6d"        "Enah"         
## [173] "Cd36"          "Diras2"        "Cd9"           "Alas2"        
## [177] "Cxcl14"        "H2-Q7"         "Sag"           "D630023F18Rik"
## [181] "Gnaz"          "Ifi30"         "H2-D1"         "C3"           
## [185] "Tspo"          "Dusp2"         "Ccl9"          "Ifitm2"       
## [189] "Gm43409"       "Ebf1"          "Fcrla"         "Spaar"        
## [193] "Lrr1"          "Mcemp1"        "Ptcra"         "Gbp2"         
## [197] "Ifi213"        "Lgals1"        "Ifi209"        "Vcam1"        
## [201] "Syde1"         "Fpr2"          "Myl4"          "Nes"          
## [205] "Slc4a4"        "Clec12a"       "Mt1"           "Bcl2a1d"      
## [209] "mt-Nd1"        "Jsrp1"         "Tpx2"          "Bend4"        
## [213] "Arntl2"        "4930414N06Rik" "Cenpe"         "Ifi206"       
## [217] "C3ar1"         "Cd40"          "Slfn1"         "Parp14"       
## [221] "Marcksl1"      "Sox2"          "Scg3"          "Gm4951"       
## [225] "Herpud1"       "Gem"           "B630019A10Rik" "Hspa1a"       
## [229] "Ebf4"          "Lrrc31"        "Adgrg5"        "Cd14"         
## [233] "Neb"           "Cxcr2"         "Ifitm1"        "Edn1"         
## [237] "Ncl"           "Hspa8"         "Sdc4"          "Nfkbia"       
## [241] "Apob"          "Id2"           "Junb"          "Pcgf2"        
## [245] "Igkc"          "Nfkbiz"        "Irf7"          "Ccna2"        
## [249] "Cmpk2"         "Fgl2"          "Ctsc"          "Lgals3bp"     
## [253] "Aoah"          "Chchd10"       "AA467197"      "Cst7"         
## [257] "B230312C02Rik" "Slfn2"         "Kif4"          "Birc5"        
## [261] "mt-Co3"        "Ly6i"          "S100a8"        "Zfp300"       
## [265] "Ifi208"        "Mzb1"          "P2ry12"        "Tnip3"        
## [269] "Ccl6"          "Gpr18"         "Ifit1bl1"      "Oas3"         
## [273] "Sox11"         "Dnaja1"        "3830403N18Rik" "mt-Atp6"      
## [277] "Il4i1"         "Gpr65"         "Cpm"           "Ecrg4"        
## [281] "Ehf"           "Gnb4"          "She"           "Neurod6"      
## [285] "Gm44684"       "Retn"          "Hp"            "Ppfia2"       
## [289] "Olfr1393"      "Mfap4"         "Gm4755"        "Ccdc88c"      
## [293] "Dlg5"          "Mov10l1"       "Sst"           "Cd200"        
## [297] "Rit2"          "Ccdc87"        "Gm1141"        "Col5a1"

PCA

# exclude MT genes  and more 
#   add sex-related Xist/Tsix

DIG <- grep("^Tra|^Trb|^Trg|^Trd|^Tcr|^Igm|^Igh|^Igk|^Igl|Jchain|Mzb1|Vpreb|Lars2|Jun|Fos|Egr1|^Hsp|^Rps|^Rpl|Hbb-|Hba-|^Dnaj|^AC|^AI|^AA|^AW|^AY|^BC|^Gm|^Hist|Rik$|-ps|Xist|Tsix|^Ifi|^Isg|^Mcm",
            rownames(GEX.seur),value = T)
CC_gene <- Hmisc::capitalize(tolower(as.vector(unlist(cc.genes.updated.2019))))
MT_gene <- grep("^mt-",rownames(GEX.seur),value = T)

VariableFeatures(GEX.seur) <- setdiff(VariableFeatures(object = GEX.seur),
                                                c(MT_gene,
                                                  DIG, 
                                                  CC_gene) )

GEX.seur <- RunPCA(GEX.seur, features = VariableFeatures(GEX.seur))
## PC_ 1 
## Positive:  Ccl12, Ccl2, Bst2, Slfn2, Ms4a6c, Gpr84, Lgals3bp, Ctsc, Fth1, Fcgr4 
##     Tspo, Srgn, Rtp4, Nme2, Calr, Fcgr1, Slfn5, Ms4a6d, Sdf2l1, Naaa 
##     Cxcl16, Manf, Ly6e, C5ar1, Ctsl, Cybb, Trim30a, Cd72, Oasl2, Clec2d 
## Negative:  P2ry12, Gpr34, Ctsd, Hpgd, Cd9, Ltc4s, Arhgap5, Sox4, Slc2a5, Pmp22 
##     Crybb1, Havcr2, Pmepa1, Slc40a1, Klhl24, Sgk1, Lst1, Ddit4, Il7r, Fam102b 
##     Thrsp, Ramp1, Slc12a2, Ttc28, Zbtb20, Tnfrsf17, Cbfa2t3, Ccl6, Crebrf, Tjp1 
## PC_ 2 
## Positive:  Ch25h, Cybb, Ccl3, Ccl4, Aoah, Cp, Ccl5, Fpr1, Saa3, Cd72 
##     Cxcl10, Cd69, Ehd1, Acod1, Clec4a1, Marco, Vcam1, Ccl7, Lilrb4a, Cxcl16 
##     Rsad2, Oasl1, Gpr18, Gadd45b, Fpr2, Tnfaip2, Sod2, Mcemp1, Bcl2a1d, Tlr2 
## Negative:  Spint1, Il4ra, Gas6, Ctsl, Ccr5, Ly6e, Gbp7, Slfn2, Atp6v0a2, Gadd45g 
##     Naaa, Fcgr2b, Sdf2l1, Fcgr1, Cd244a, Calr, Cdc42ep3, Ccl12, Lifr, Cct6a 
##     Trim30a, Nav3, Tuba1a, Manf, Bhlhe41, Srgn, Nme2, Socs3, Igfbp6, Pdia3 
## PC_ 3 
## Positive:  Rsad2, Iigp1, Usp18, Oasl2, Slfn5, Oasl1, Cxcl10, Cmpk2, Irf7, Fgl2 
##     Herc6, Tnfsf10, Parp14, Gbp2, Samd9l, Phf11d, Slfn9, Rnf213, Gbp5, Tlr3 
##     Ddx58, Samhd1, Slfn8, Stat1, Phf11b, Tgtp2, Igtp, Irgm1, Stat2, Sp100 
## Negative:  Npm1, Msr1, Ehd1, Eef1g, Ncl, Ran, Lpl, Anp32b, Gnl3, Tmsb4x 
##     Lilrb4a, Nop56, Eif5a, Cpd, Tlr2, Ftl1, C5ar1, Gpr84, Nme2, Ranbp1 
##     Sod2, Nop58, Pilrb1, Nme1, Tnf, Bola2, Il1b, Ctsz, Fpr1, Eif2s2 
## PC_ 4 
## Positive:  Fcer1g, Tmsb4x, Fth1, Ly86, Ftl1, B2m, Apoe, Lyz2, Cd52, Tspo 
##     H2-D1, Ctsl, Cd63, Bst2, Nme2, Ssr4, Tmem176a, Actb, Prdx5, Tmem176b 
##     H2-K1, Psme1, Atp5g1, Zbp1, Psme2, Cst7, Spp1, Gapdh, Ly6a, Cox7b 
## Negative:  Macf1, Kcnq1ot1, Dst, Zeb2, Filip1l, Arhgap5, Tnf, Fam102b, Hivep3, Ccl4 
##     Pou2f2, Fchsd2, Ranbp2, Trip11, Slc15a3, Ankrd11, Il1b, Zfp36l1, Nav3, Tnfaip2 
##     Tlr2, Atp2a2, C3ar1, Baz1b, Ncl, Runx1, Tjp1, Cep350, Btg2, Sik2 
## PC_ 5 
## Positive:  Cebpd, Il1b, Ccl12, Herpud1, Tnf, Ly6e, Ch25h, Fcgr1, Cebpb, Socs3 
##     C5ar1, Manf, Pdia6, Gpr84, Ms4a6d, Tmem176b, Slc2a5, Tnfaip2, Ms4a6b, Srgn 
##     Ms4a6c, Ccl2, Il1a, Ehd1, Sdf2l1, Saa3, Tuba1a, Gpr18, Fam102b, Nfkbia 
## Negative:  Apoe, Ctsb, Chd4, Ptma, B2m, Ftl1, Lpl, Gatm, Pkm, Cd52 
##     Stmn1, Lgals3, Itgax, Npnt, Scd2, Axl, Cstb, Mt1, Abcg1, Emp3 
##     Igfbp4, Lyz2, H2-D1, Atrx, Macf1, Actb, Npm1, Aplp2, Nap1l1, Lmnb1
#grep("tdTomato|Spp1|Cre",VariableFeatures(GEX.seur),value = T)
length(VariableFeatures(GEX.seur))
## [1] 1437
head(VariableFeatures(GEX.seur),500)
##   [1] "Saa3"       "Cxcl10"     "Ccl4"       "Ccl3"       "Ccl5"      
##   [6] "Marco"      "Cxcl13"     "Apoe"       "Ch25h"      "Lcn2"      
##  [11] "Ccl2"       "Ccl7"       "Ccl12"      "Spp1"       "Rsad2"     
##  [16] "Il1b"       "Cd69"       "Cxcl9"      "Acod1"      "Slc13a3"   
##  [21] "Cybb"       "Iigp1"      "Il12b"      "Cd74"       "Lpl"       
##  [26] "Lyz2"       "Tnf"        "Cp"         "Egr3"       "Cd72"      
##  [31] "Fth1"       "Fpr1"       "Ly6k"       "Aldh1a3"    "Gpr84"     
##  [36] "Wfdc21"     "Olfr889"    "H2-Ab1"     "Lgals3"     "H2-Aa"     
##  [41] "Wfdc17"     "Ccr7"       "Cxcl2"      "Fn1"        "Ms4a6c"    
##  [46] "Atf3"       "Loxl2"      "Tnfaip2"    "Prc1"       "Dennd5b"   
##  [51] "Wdcp"       "H2-Eb1"     "Cd52"       "Pbld1"      "Kcnk4"     
##  [56] "Gdf15"      "Oasl1"      "Gadd45b"    "Ptges"      "Cxcl16"    
##  [61] "Egr2"       "Ccr1"       "Ptgds"      "Fabp5"      "Apoc1"     
##  [66] "Cd79a"      "Ms4a4c"     "Tlr2"       "Dqx1"       "Il10"      
##  [71] "Esco2"      "Saa1"       "E2f7"       "Slfn5"      "Cpd"       
##  [76] "Nr4a1"      "Lockd"      "Bst2"       "Lilrb4a"    "Kcnq1ot1"  
##  [81] "Pclaf"      "Csf1"       "Cd63"       "Pbk"        "Ier3"      
##  [86] "Cd83"       "Mir155hg"   "Sod2"       "Clec4a1"    "Lhfpl4"    
##  [91] "Ptprcap"    "Rhox5"      "Cdkn1a"     "Msr1"       "Dnase1l3"  
##  [96] "Ifnb1"      "Plk1"       "C4b"        "Stmn1"      "Srgn"      
## [101] "Spib"       "Ly6a"       "H2-K1"      "Vgf"        "Slc1a2"    
## [106] "Dennd3"     "Vnn3"       "Oasl2"      "Plac8"      "Usp18"     
## [111] "Mmp12"      "Ptn"        "Fcgr4"      "Ms4a6d"     "Enah"      
## [116] "Cd36"       "Diras2"     "Cd9"        "Alas2"      "Cxcl14"    
## [121] "H2-Q7"      "Sag"        "Gnaz"       "H2-D1"      "C3"        
## [126] "Tspo"       "Dusp2"      "Ccl9"       "Ebf1"       "Fcrla"     
## [131] "Spaar"      "Lrr1"       "Mcemp1"     "Ptcra"      "Gbp2"      
## [136] "Lgals1"     "Vcam1"      "Syde1"      "Fpr2"       "Myl4"      
## [141] "Nes"        "Slc4a4"     "Clec12a"    "Mt1"        "Bcl2a1d"   
## [146] "Jsrp1"      "Bend4"      "Arntl2"     "C3ar1"      "Cd40"      
## [151] "Slfn1"      "Parp14"     "Marcksl1"   "Sox2"       "Scg3"      
## [156] "Herpud1"    "Gem"        "Ebf4"       "Lrrc31"     "Adgrg5"    
## [161] "Cd14"       "Neb"        "Cxcr2"      "Edn1"       "Ncl"       
## [166] "Sdc4"       "Nfkbia"     "Apob"       "Id2"        "Pcgf2"     
## [171] "Nfkbiz"     "Irf7"       "Ccna2"      "Cmpk2"      "Fgl2"      
## [176] "Ctsc"       "Lgals3bp"   "Aoah"       "Chchd10"    "Cst7"      
## [181] "Slfn2"      "Kif4"       "Ly6i"       "S100a8"     "Zfp300"    
## [186] "P2ry12"     "Tnip3"      "Ccl6"       "Gpr18"      "Oas3"      
## [191] "Sox11"      "Il4i1"      "Gpr65"      "Cpm"        "Ecrg4"     
## [196] "Ehf"        "Gnb4"       "She"        "Neurod6"    "Retn"      
## [201] "Hp"         "Ppfia2"     "Olfr1393"   "Mfap4"      "Ccdc88c"   
## [206] "Dlg5"       "Mov10l1"    "Sst"        "Cd200"      "Rit2"      
## [211] "Ccdc87"     "Col5a1"     "Sel1l3"     "Mmd2"       "Plbd1"     
## [216] "Spock3"     "Alox8"      "Mdga2"      "Chrm3"      "Synpr"     
## [221] "Ubash3a"    "Ank"        "Milr1"      "Ehd1"       "Il2rg"     
## [226] "Rnf213"     "Zbtb20"     "Adra2a"     "Herc6"      "Mrc1"      
## [231] "Ctla2a"     "Clec2d"     "Gpr34"      "Manf"       "Samd9l"    
## [236] "Trim17"     "Prdx5"      "Il1a"       "Xdh"        "Lox"       
## [241] "Akap12"     "Il1rn"      "Ly86"       "Tpr"        "Phf11b"    
## [246] "Fam20c"     "Syne1"      "Usp26"      "Enpp2"      "Cobll1"    
## [251] "Ocstamp"    "Robo2"      "Calr"       "Ctsl"       "Klf2"      
## [256] "Npm1"       "Pcp4l1"     "Ptprz1"     "Smc2"       "Phf11d"    
## [261] "Gbp5"       "Cnn3"       "Spdya"      "Sdf2l1"     "Gbp7"      
## [266] "Sox4"       "Pilra"      "Gbp4"       "Apba2"      "Ahnak2"    
## [271] "Gap43"      "Insl6"      "Slfn9"      "Olfr1423"   "Rtp4"      
## [276] "Igsf6"      "Plcb1"      "Dpp10"      "Tmeff1"     "A3galt2"   
## [281] "Tcea3"      "Kif17"      "Srrm4"      "Cttnbp2"    "Napsa"     
## [286] "Tulp2"      "Zfp941"     "Olfr921"    "Scml4"      "Gipc3"     
## [291] "Cfap54"     "Rflnb"      "Olfr113"    "Actn3"      "Myof"      
## [296] "Rab33a"     "Gpc4"       "Brs3"       "Fndc3c1"    "Sgk1"      
## [301] "Ptger4"     "Ms4a7"      "Mif"        "Gnao1"      "C5ar1"     
## [306] "Bicd1"      "Anks1b"     "Eif2ak2"    "Stat1"      "Rad50"     
## [311] "Morc4"      "Aox2"       "Pm20d1"     "Nfasc"      "Zfp385b"   
## [316] "Tshz2"      "Bhlhe22"    "Igsf3"      "Synpo2"     "Dkk2"      
## [321] "Ptprf"      "Pcdh7"      "Lmntd2"     "Pou2af1"    "Ucn2"      
## [326] "Gli1"       "Slc13a5"    "Ccdc92b"    "Slfn4"      "Pcdh8"     
## [331] "Pou4f1"     "Csdc2"      "A4galt"     "Lipo2"      "Nyx"       
## [336] "Lpar4"      "Lst1"       "Il7r"       "Cfb"        "Trim30a"   
## [341] "Brca1"      "Plp1"       "Ccdc163"    "Meis2"      "Poln"      
## [346] "H2-Q6"      "Nme2"       "Adgre1"     "Ier2"       "Naaa"      
## [351] "Btg2"       "Pdlim7"     "Actr3b"     "Spock1"     "Bcl2a1a"   
## [356] "Ankrd12"    "Ms4a6b"     "Cep290"     "Rnf217"     "Strip2"    
## [361] "Nrgn"       "Clu"        "Smc1b"      "Ccdc34"     "Ikbke"     
## [366] "Cox6a2"     "Ctsb"       "Taf1d"      "Ftl1"       "Itgal"     
## [371] "Nhs"        "Ctsd"       "Eif5b"      "Cct6b"      "Ccl8"      
## [376] "Crybb1"     "Upk1b"      "Plaur"      "Icam4"      "Iqch"      
## [381] "Tceal3"     "Pmp22"      "Slc31a2"    "Ccdc88a"    "Lilr4b"    
## [386] "Adamts1"    "Arpp21"     "Aspm"       "Clec2i"     "Mpp4"      
## [391] "Mpz"        "Kcnk2"      "Bmp2"       "Schip1"     "Fam110b"   
## [396] "Galnt9"     "Ffar1"      "Kcnc3"      "Ctrl"       "Cdon"      
## [401] "Tmem25"     "Rbms3"      "Grip1"      "Slc25a21"   "Arg2"      
## [406] "Gpr132"     "Spata31d1a" "Adra1a"     "Nefm"       "Pcdhb4"    
## [411] "Pcdhb7"     "Rhod"       "Plk4"       "Il4ra"      "Gapdh"     
## [416] "Axl"        "Spats2l"    "Ly6e"       "Ldlr"       "Zbp1"      
## [421] "Map7d3"     "Pilrb2"     "Rcan1"      "Eif3a"      "Socs3"     
## [426] "Cxcr4"      "Nrxn1"      "Socs1"      "Klf4"       "Sgo2a"     
## [431] "Brip1"      "Pilrb1"     "Sgo1"       "Sparcl1"    "Ddx58"     
## [436] "Ccnb1"      "Fcgr2b"     "Dmrtb1"     "Mns1"       "Mt2"       
## [441] "Sp100"      "Phlda3"     "Tor3a"      "Smc1a"      "Mndal"     
## [446] "Cspg4"      "Hcar2"      "Knl1"       "Ankrd35"    "Spata7"    
## [451] "Zfp811"     "Agbl2"      "Rragd"      "Xaf1"       "Ddx21"     
## [456] "Etaa1os"    "Runx1t1"    "Ccnd1"      "Fcgr1"      "Fmr1nb"    
## [461] "Ankrd26"    "Akap9"      "Eea1"       "Postn"      "Cadm1"     
## [466] "Socs2"      "Ezh2"       "Samhd1"     "Slamf8"     "Emp3"      
## [471] "Golgb1"     "Icam1"      "Ccdc63"     "Olfr1328"   "Crlf2"     
## [476] "Tmem52"     "Efna5"      "Pclo"       "Cldn7"      "Nefl"      
## [481] "Pcp4"       "Dusp1"      "Knstrn"     "Anp32b"     "Tmem176a"  
## [486] "Ranbp1"     "Morn5"      "Kcnip4"     "Smco3"      "Adgre5"    
## [491] "Tgtp1"      "Rab15"      "Prrt1"      "Daam2"      "Kcnip2"    
## [496] "Cldn34b2"   "Igfbp5"     "Itga8"      "Ly75"       "Kcnh7"
DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "orig.ident") +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "orig.ident")

DimPlot(GEX.seur, reduction = "pca",dims = 1:2, group.by = "FB.info", cols = color.FB2) +
  DimPlot(GEX.seur, reduction = "pca",dims = 3:4, group.by = "FB.info", cols = color.FB2)

DimHeatmap(GEX.seur, dims = 1:12, cells = 1500, balanced = TRUE,ncol = 4)

(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(desc(PC_1) ) %>% head(40)
##                PC_1         PC_2          PC_3         PC_4         PC_5
## Ccl12    0.14457677 -0.058118545 -0.0018318181  0.007582612  0.082396421
## Ccl2     0.11499884  0.063417638  0.0017777238 -0.036221048  0.051076874
## Bst2     0.11287944 -0.021007628  0.0277652735  0.076082452  0.010035533
## Slfn2    0.10685377 -0.071445499  0.0401899236  0.003519736  0.014894258
## Ms4a6c   0.10574457  0.038831761 -0.0014509955  0.014226352  0.052684192
## Gpr84    0.10539918  0.024039318 -0.0540462848 -0.058538073  0.063324958
## Lgals3bp 0.10417045 -0.049830290  0.0362641497  0.039165144 -0.011137749
## Ctsc     0.10095294 -0.025918417 -0.0200639846  0.007531145  0.042178753
## Fth1     0.09876637  0.050566857 -0.0450666237  0.122084274  0.018421379
## Fcgr4    0.09828840 -0.001948747  0.0240950665  0.048042450  0.021747200
## Tspo     0.09827368 -0.048911736  0.0275754640  0.092588784 -0.002869758
## Srgn     0.09716155 -0.054021370 -0.0152596736  0.030543617  0.052736727
## Rtp4     0.09018160 -0.048992222  0.0697324906  0.034217045  0.014267289
## Nme2     0.08832746 -0.053470280 -0.0531938637  0.073291733 -0.026093787
## Calr     0.08651723 -0.059477349 -0.0139172136 -0.016849079  0.040639184
## Fcgr1    0.08586060 -0.061155353  0.0478785674  0.017247831  0.074430496
## Slfn5    0.08555578 -0.006093494  0.1361551966  0.009873710 -0.026867442
## Ms4a6d   0.08544524  0.011354597 -0.0038894255  0.026874812  0.061624398
## Sdf2l1   0.08221451 -0.062916160 -0.0151158393  0.043751016  0.048960745
## Naaa     0.08077687 -0.066256678 -0.0060345557  0.023355379  0.009547661
## Cxcl16   0.08052397  0.092534600 -0.0473348881 -0.018873252  0.008203213
## Manf     0.08045464 -0.055168054 -0.0194439437  0.013993754  0.066358700
## Ly6e     0.08004240 -0.076720221  0.0467876609  0.047706314  0.074941955
## C5ar1    0.07987303 -0.052252887 -0.0549849198 -0.040784930  0.067407259
## Ctsl     0.07982985 -0.080596382 -0.0261996879  0.088109352 -0.005454459
## Cybb     0.07939603  0.167123491 -0.0126246767 -0.005137725  0.014906403
## Trim30a  0.07689381 -0.056496199  0.0823487809 -0.000542407  0.011297798
## Cd72     0.07648303  0.125824294 -0.0388548602  0.018715671  0.010916926
## Oasl2    0.07646800  0.013691289  0.1524426284  0.034955530 -0.028014516
## Clec2d   0.07581392  0.041383796  0.0422711132 -0.013785900  0.027425998
## Pdia3    0.07481310 -0.052545826 -0.0033332953 -0.023359192  0.013708430
## Ssr4     0.07423175 -0.025851064 -0.0454374940  0.071868577 -0.006559799
## Msr1     0.07401689  0.050445457 -0.0657520967 -0.047195511  0.019780185
## Cebpb    0.07222678 -0.048136119 -0.0136726779 -0.020023584  0.074376991
## Ccdc86   0.07117643 -0.033635235 -0.0006251396 -0.014705320  0.003487167
## Zbp1     0.07115050  0.005481147  0.0840959043  0.054776163 -0.014157260
## Pdia6    0.07020103 -0.047932436 -0.0154576097 -0.005955612  0.064774634
## Rnf213   0.06963596 -0.006098082  0.1007877602 -0.015361398 -0.028253033
## Usp18    0.06959384  0.022581800  0.1554107871  0.002425291 -0.011111606
## Ch25h    0.06908844  0.202951394 -0.0223728688 -0.056313752  0.074726974
##                   PC_6          PC_7         PC_8
## Ccl12     0.0136698459  0.0046974571  0.064177717
## Ccl2      0.0379487022  0.0632854552  0.039182200
## Bst2      0.0070299225 -0.0188225149 -0.042249980
## Slfn2    -0.0063884225  0.0095451345  0.049317330
## Ms4a6c   -0.0317392773 -0.0650964964 -0.012636281
## Gpr84     0.0692366858 -0.0025107600  0.076305163
## Lgals3bp -0.0072603718 -0.0738942966  0.002982274
## Ctsc      0.0178042744 -0.0620925892  0.009439878
## Fth1      0.0251511785 -0.0631068145 -0.035943539
## Fcgr4    -0.0150690086 -0.0597129164  0.003799063
## Tspo     -0.0028361490 -0.0070424635  0.011198939
## Srgn      0.0745857920 -0.0513106858 -0.024158239
## Rtp4     -0.0044163405 -0.0002538197  0.004158800
## Nme2     -0.0025394374  0.0497368308 -0.020111354
## Calr      0.0098756177 -0.0692140482 -0.072194278
## Fcgr1     0.0004209819 -0.0124601044  0.018466400
## Slfn5     0.0235324494 -0.0090961311  0.022404268
## Ms4a6d   -0.0252149058 -0.0669575592 -0.041855782
## Sdf2l1    0.0126317948 -0.0564524766 -0.092176552
## Naaa     -0.0297093839 -0.0199996907  0.042660378
## Cxcl16    0.0630800838 -0.0556722309  0.015232151
## Manf      0.0100804161 -0.0376948029 -0.079910284
## Ly6e      0.0327829853 -0.0632113534 -0.010175082
## C5ar1     0.0681588607 -0.0003147902  0.032370918
## Ctsl      0.0060387839 -0.0081628936  0.058153278
## Cybb     -0.0623195004 -0.0477167605 -0.033990247
## Trim30a  -0.0015685118 -0.0283835209  0.008363304
## Cd72     -0.0902896994 -0.0763147539 -0.021445219
## Oasl2     0.0073778915 -0.0054387363 -0.003204502
## Clec2d   -0.0069487766 -0.0038654131  0.053390801
## Pdia3    -0.0028536264 -0.0852722311 -0.073132538
## Ssr4      0.0112912699  0.0080203216 -0.031416300
## Msr1     -0.0470665789 -0.0197675261 -0.036737132
## Cebpb     0.0448126464  0.0354570713  0.035654950
## Ccdc86    0.0193367799  0.0253798051  0.001065169
## Zbp1      0.0025978752 -0.0067053024 -0.023801228
## Pdia6     0.0151204901 -0.0581104673 -0.083242720
## Rnf213    0.0076861271 -0.0184715440 -0.003721215
## Usp18     0.0121419387  0.0662070870 -0.042097205
## Ch25h    -0.0944706090  0.0476184573 -0.082918648
(GEX.seur@reductions$pca@feature.loadings %>% as.data.frame)[,1:8] %>% arrange(PC_1 ) %>% head(40)
##                 PC_1         PC_2         PC_3         PC_4          PC_5
## P2ry12   -0.13273384 -0.005629467  0.065485216 -0.058335671  0.0151140854
## Gpr34    -0.12544252  0.015181896  0.038844037 -0.037886241  0.0028173495
## Ctsd     -0.08913450 -0.002412970  0.031616554  0.017722474  0.0255570136
## Hpgd     -0.07686498  0.008751071  0.022677487 -0.015074420 -0.0211753874
## Cd9      -0.06941297  0.079331428 -0.025401577 -0.009565905 -0.0391111457
## Ltc4s    -0.06347304 -0.037181944  0.012272317  0.015825692  0.0042799467
## Arhgap5  -0.06343962 -0.021620906  0.042195911 -0.087940154  0.0004442022
## Sox4     -0.06032027  0.019189153  0.002317105 -0.038065717 -0.0184256103
## Slc2a5   -0.05778217 -0.017066611  0.034648291 -0.045898963  0.0554507374
## Pmp22    -0.05631634  0.010509722  0.001787813 -0.056270124  0.0295174248
## Crybb1   -0.05363518  0.003096415  0.019754398  0.031836780 -0.0555371069
## Havcr2   -0.05241728  0.025417142  0.020895802 -0.048990107 -0.0103852553
## Pmepa1   -0.05241534 -0.010455675  0.038619220 -0.038216325  0.0121668260
## Slc40a1  -0.04987126  0.006024600  0.030980473 -0.058748445  0.0051023656
## Klhl24   -0.04734085  0.027338254  0.012989055 -0.024154349 -0.0196886328
## Sgk1     -0.04631263  0.023857655  0.004370984 -0.037970720 -0.0226033893
## Lst1     -0.04405431  0.056891358  0.006394751  0.039791783 -0.0229769178
## Ddit4    -0.04263706  0.004365984  0.021918562 -0.033052296  0.0316676731
## Il7r     -0.04238146  0.010493410  0.037419550 -0.030324366  0.0310134946
## Fam102b  -0.04081603 -0.013585604  0.016690231 -0.084004907  0.0468995186
## Thrsp    -0.04073298  0.015527898  0.019420394 -0.015585311  0.0432943542
## Ramp1    -0.04048046  0.031555175 -0.006513382  0.017861894 -0.0680257710
## Slc12a2  -0.04007155  0.017772048  0.005035451 -0.037857276 -0.0578854071
## Ttc28    -0.03850200  0.012653676  0.018642935 -0.041098095 -0.0195166438
## Zbtb20   -0.03793839  0.009000660  0.031279739 -0.042489195  0.0074691960
## Tnfrsf17 -0.03626810  0.020521373  0.010148565 -0.016869586  0.0140088553
## Cbfa2t3  -0.03488123  0.003437202  0.011728702 -0.021272566 -0.0170343673
## Ccl6     -0.03406132  0.025313415 -0.003239522  0.018979392 -0.0433961278
## Crebrf   -0.03207459  0.018683094  0.024630248 -0.016528276 -0.0206936351
## Tjp1     -0.03200141 -0.004899833  0.022942055 -0.063871496 -0.0181313773
## Rapsn    -0.03076310  0.023023012  0.005218324  0.010725474 -0.0160285699
## Upk1b    -0.03022652  0.005547926  0.015543692 -0.037821669  0.0458513738
## Arid5b   -0.02982933  0.019527703  0.030022740 -0.018714203 -0.0188086965
## Tent5c   -0.02871037  0.012815758  0.008879275 -0.032698148 -0.0103643085
## S1pr1    -0.02846519  0.023880914 -0.004418825 -0.012039410 -0.0330202841
## Foxp1    -0.02803479  0.003406680  0.015871555 -0.015308127 -0.0132311851
## Csmd3    -0.02775643 -0.017580578  0.017069868 -0.061550180  0.0341682724
## Fry      -0.02731603  0.005309651  0.019753401 -0.032771511 -0.0228183670
## Hacd4    -0.02722622  0.019081804  0.017959105  0.018999979 -0.0224948846
## Ppm1l    -0.02714267  0.003988630  0.015116721 -0.031448084  0.0114350924
##                  PC_6         PC_7          PC_8
## P2ry12   -0.026156982  0.018240462 -0.0300755193
## Gpr34     0.026532943 -0.006372836 -0.0651009026
## Ctsd      0.190974686 -0.104487648 -0.1260501454
## Hpgd     -0.031284412  0.049200663 -0.0083889914
## Cd9       0.151841316 -0.024208949 -0.1155900492
## Ltc4s    -0.012277965  0.084374302  0.0178319841
## Arhgap5  -0.019241640 -0.059059787 -0.0760211461
## Sox4      0.011739737  0.020232577 -0.0432704868
## Slc2a5   -0.029081599  0.013571293 -0.0739490849
## Pmp22     0.110941775 -0.025663507 -0.1594533191
## Crybb1   -0.049237065  0.070428344  0.0243721919
## Havcr2    0.049704513 -0.063001520 -0.0621450743
## Pmepa1   -0.004498927  0.001010207 -0.0675075488
## Slc40a1  -0.007718383 -0.036141032 -0.0584239484
## Klhl24   -0.002514213 -0.030520337 -0.0049406889
## Sgk1      0.039579319 -0.002879153 -0.0440935219
## Lst1     -0.044907219 -0.001087295  0.0260939418
## Ddit4    -0.024443147 -0.002190533 -0.0586518156
## Il7r     -0.009335944 -0.031865130 -0.0340257081
## Fam102b   0.008957638 -0.013137921 -0.0930676820
## Thrsp     0.024279140 -0.020148595 -0.0731248365
## Ramp1    -0.030205746  0.046670420  0.0623262694
## Slc12a2   0.011754490 -0.013417529 -0.0304803012
## Ttc28     0.009818731 -0.029149891 -0.0382097248
## Zbtb20   -0.001513918 -0.041651298 -0.0199665265
## Tnfrsf17 -0.001118765  0.021941105 -0.0322695726
## Cbfa2t3  -0.018356741  0.013471296 -0.0007691025
## Ccl6      0.051207150 -0.005939162 -0.0448838530
## Crebrf   -0.008592314 -0.021871757 -0.0151609124
## Tjp1     -0.022404650 -0.022099382 -0.0401320487
## Rapsn     0.024991789 -0.011592548 -0.0407302242
## Upk1b     0.009742560  0.018709700 -0.0610980231
## Arid5b    0.010423645 -0.030036734 -0.0314682264
## Tent5c    0.060658674 -0.037930916 -0.0659988556
## S1pr1     0.063573154 -0.029272221 -0.0560909089
## Foxp1    -0.024303502 -0.003825656 -0.0012504121
## Csmd3    -0.009686347  0.021762739 -0.0748941670
## Fry      -0.032936399 -0.006355581  0.0046956379
## Hacd4    -0.002810830 -0.007709384  0.0328864595
## Ppm1l    -0.009549849 -0.001740995 -0.0240216106
decide PCs to use
ElbowPlot(GEX.seur,ndims = 50)

PCs <- 1:20
GEX.seur <- FindNeighbors(GEX.seur, dims = PCs, k.param = 20)
## Computing nearest neighbor graph
## Computing SNN
GEX.seur <- FindClusters(GEX.seur, method = 'igraph' ,resolution = 1.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 15536
## Number of edges: 394924
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6553
## Number of communities: 15
## Elapsed time: 2 seconds

Run UMAP/tSNE

GEX.seur <- RunTSNE(GEX.seur, dims=PCs, complexity = 100)
GEX.seur <- RunUMAP(GEX.seur, dims=PCs, n.neighbors = 20, seed.use = 998)
## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session
## 22:36:25 UMAP embedding parameters a = 0.9922 b = 1.112
## 22:36:25 Read 15536 rows and found 20 numeric columns
## 22:36:25 Using Annoy for neighbor search, n_neighbors = 20
## 22:36:25 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 22:36:27 Writing NN index file to temp file C:\Users\Shaorui\AppData\Local\Temp\RtmpI5QH6X\filed2d056b0165
## 22:36:27 Searching Annoy index using 1 thread, search_k = 2000
## 22:36:31 Annoy recall = 100%
## 22:36:32 Commencing smooth kNN distance calibration using 1 thread
## 22:36:33 Initializing from normalized Laplacian + noise
## 22:36:33 Commencing optimization for 200 epochs, with 451104 positive edges
## 22:36:49 Optimization finished
DimPlot(GEX.seur, reduction = "tsne", label = T) + DimPlot(GEX.seur, reduction = "umap", label = T)

FeaturePlot(GEX.seur, reduction = "umap", features = c("nFeature_RNA","nCount_RNA","percent.mt","percent.rb"))

GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
                          levels = levels(GEX.seur$FB.info)[c(9,8,6,7,5:3)])
DimPlot(GEX.seur, label = F, pt.size = 0.05, repel = F, reduction = 'umap', group.by = "FB.new", split.by = "FB.new", 
        ncol = 4, cols = color.FB2[c(7,6,4,5,3,2,1)])

DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters") +
  DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.info", cols = color.FB2)

check some markers

markers.mig <- c("Ptprc",#"Cd3d","Cd3e","Cd19",
                 "Cd74","Lyz2","Ccl4",
                 "Aif1","P2ry12","C1qa","Spp1",
                 "Top2a","Pcna","Mki67","Mcm6",
                 "Cx3cr1","Il4ra","Il13ra1","Fabp5",
                 "Hmox1","Ms4a7","Pf4","Tubb3",
                 "tdTomato-head","tdTomato-mid","tdTomato-tail","CreERT2"
                 )
FeaturePlot(GEX.seur, 
            features = markers.mig,
            ncol = 4)

DotPlot(GEX.seur, features = markers.mig, group.by = "FB.info",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

VlnPlot(GEX.seur, features = c("tdTomato-head","tdTomato-mid","tdTomato-tail","CreERT2"),
        group.by = "FB.info", ncol = 2)

DotPlot(GEX.seur, features = c(markers.mig,"Ifit1","Ifit2"), group.by = "seurat_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

sort

GEX.seur$sort_clusters <- factor(as.character(GEX.seur$seurat_clusters),
                                 levels = c(2,1,8,7,9,
                                            14,
                                            12,3,5,0,4,6,11,13,10))
cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$FB.new,
      clusters=GEX.seur$sort_clusters),
                   main = "Cell Count",
                   gaps_row = c(4),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$FB.new,
                                clusters=GEX.seur$sort_clusters)),
                   main = "Cell Ratio",
                   gaps_row = c(4),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)

markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "sort_clusters"

#GEX.markers.pre <- FindAllMarkers(GEX.seur, only.pos = TRUE, 
#                                  min.pct = 0.05,
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.25)
GEX.markers.pre <- read.table("sc10x_LYN.marker.subset2_sort1007.csv", header = TRUE, sep = ",")
GEX.markers.pre %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 120 x 7
## # Groups:   cluster [15]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  
##        <dbl>      <dbl> <dbl> <dbl>     <dbl>   <int> <chr> 
##  1 0              0.887 0.989 0.922 0               2 Gpr34 
##  2 0              0.858 1     0.973 0               2 P2ry12
##  3 5.16e-226      0.604 0.981 0.963 9.24e-222       2 Rgs10 
##  4 1.17e-204      0.652 0.921 0.855 2.09e-200       2 Rnase4
##  5 4.07e-203      0.707 0.923 0.804 7.29e-199       2 Cd9   
##  6 1.13e-192      0.724 0.796 0.665 2.03e-188       2 Hpgd  
##  7 1.56e-143      0.603 0.586 0.48  2.80e-139       2 Ang   
##  8 4.46e-130      0.613 0.583 0.477 7.98e-126       2 Lst1  
##  9 4.83e-170      0.923 0.962 0.847 8.64e-166       7 H2-D1 
## 10 5.05e-108      0.582 0.996 0.926 9.04e-104       7 Gpr34 
## # ... with 110 more rows
GEX.markers.sort <- GEX.markers.pre
GEX.markers.sort$cluster <- factor(as.character(GEX.markers.sort$cluster),
                          levels = levels(GEX.seur$sort_clusters))


markers.pre_t48 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1) %>%
                   top_n(n = 48, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t60 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene

markers.pre_t120 <- (GEX.markers.sort %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl",GEX.markers.sort$gene,invert = T,value = T)) %>%
                   top_n(n = 120, wt = avg_log2FC) %>%
                   ungroup() %>%
  arrange(desc(avg_log2FC*pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[1:64]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[65:128]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[129:192]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[193:256]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[257:320]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[321:384]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[385:448]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[449:498]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "sort_clusters", features = rev(markers.pre_t60[499:529]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

GEX.seur$cnt <- as.character(GEX.seur$FB.info)
GEX.seur$cnt[GEX.seur$cnt %in% c("P30.LPS.TDT1","P30.LPS.TDT2")] <- "P30.LPS.TDT"

GEX.seur$cnt <- factor(GEX.seur$cnt,
                       levels = c("P30.PBS.CTR","P30.PBS.MIG","P30.PBS.TDT",
                                  "P30.LPS.CTR","P30.LPS.MIG","P30.LPS.TDT"))
color.cnt <- rev(color.FB2)[c(5:7,1:2,4)]

subset2 DEGs

DEGs.a

GEX.comp <- GEX.seur
Idents(GEX.comp) <- "cnt"
GEX.comp
## An object of class Seurat 
## 17900 features across 15536 samples within 1 assay 
## Active assay: RNA (17900 features, 1437 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.comp$cnt[1:6]
## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1 
##        P30.LPS.CTR        P30.LPS.TDT        P30.PBS.MIG        P30.LPS.TDT 
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1 
##        P30.PBS.TDT        P30.LPS.MIG 
## 6 Levels: P30.PBS.CTR P30.PBS.MIG P30.PBS.TDT P30.LPS.CTR ... P30.LPS.TDT
DEGs.a <- list()


for(mm in c("P30.PBS","P30.LPS")){
  for(nn in c("CTR","MIG")){
    DEGs.a[[paste0(mm,".TDTvs",nn)]] <- FindAllMarkers(subset(GEX.comp, subset = cnt %in% c(paste0(mm,".",nn),paste0(mm,".TDT")) & age %in% c(mm)),
                                                       assay = "RNA",
                                                       only.pos = T,
                                                       min.pct = 0.05, 
                                                       logfc.threshold = 0.1,
                                                       test.use = "MAST")
    
     
  }
}
## Calculating cluster P30.PBS.CTR
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.PBS.TDT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.PBS.MIG
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.PBS.TDT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS.CTR
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS.TDT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS.MIG
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS.TDT
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
#DEGs.a
names(DEGs.a)
## [1] "P30.PBS.TDTvsCTR" "P30.PBS.TDTvsMIG" "P30.LPS.TDTvsCTR" "P30.LPS.TDTvsMIG"
# save DEGs' table
df.DEGs.a <- data.frame()

for(nn in names(DEGs.a)){
  df.DEGs.a <- rbind(df.DEGs.a, data.frame(DEGs.a[[nn]], contrast = nn))
}

#write.csv(df.DEGs.a, "DEGs.a.subset2.csv")

DEG stat

cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
##           contrast     cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR       8
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT      29
## 3 P30.LPS.TDTvsMIG P30.LPS.TDT      51
## 4 P30.LPS.TDTvsMIG P30.LPS.MIG      17
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR       2
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT       7
## 7 P30.PBS.TDTvsMIG P30.PBS.TDT       1
## 8 P30.PBS.TDTvsMIG P30.PBS.MIG       3
pp.stat.DEG <- list()
df.DEGs.a$cluster <- factor(as.character(df.DEGs.a$cluster),
                            levels = levels(GEX.seur$cnt)[c(4:6,1:3)])
pp.stat.DEG[["a"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["a"]]

cut.padj = 0.05
cut.log2FC = 0.25#0.3   
  
cut.pct1 = 0.05

df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
##           contrast     cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR      13
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT      37
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG      20
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT      72
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR       4
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT      11
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG       4
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT       3
pp.stat.DEG[["a1"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt, name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["a1"]]

DEG plot

pp.DEGs.a <- lapply(DEGs.a, function(x){NA})
names(DEGs.a)
## [1] "P30.PBS.TDTvsCTR" "P30.PBS.TDTvsMIG" "P30.LPS.TDTvsCTR" "P30.LPS.TDTvsMIG"
DEGs.a.combine <- DEGs.a
DEGs.a.combine <- lapply(DEGs.a.combine, function(x){
  for(zz in c("P30.LPS.MIG","P30.LPS.CTR","P30.PBS.MIG","P30.PBS.CTR")){
    x[x$cluster==zz,"avg_log2FC"] <- -x[x$cluster==zz,"avg_log2FC"]
  }
  
  x
})
pp.DEGs.a$P30.LPS.TDTvsCTR <- DEGs.a.combine$P30.LPS.TDTvsCTR %>%
  mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTR"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.LPS, TDT vs CTR") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.LPS.TDTvsCTR

pp.DEGs.a$P30.LPS.TDTvsMIG <- DEGs.a.combine$P30.LPS.TDTvsMIG %>%
  mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("MIG"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.LPS, TDT vs MIG") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,1))
pp.DEGs.a$P30.LPS.TDTvsMIG

pp.DEGs.a$P30.PBS.TDTvsCTR <- DEGs.a.combine$P30.PBS.TDTvsCTR %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTR"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.PBS, TDT vs CTR") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.PBS.TDTvsCTR

pp.DEGs.a$P30.PBS.TDTvsMIG <- DEGs.a.combine$P30.PBS.TDTvsMIG %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("MIG"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.PBS, TDT vs MIG") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") #+ xlim(c(-1,2))
pp.DEGs.a$P30.PBS.TDTvsMIG

DEGs.b

GEX.comp <- GEX.seur
Idents(GEX.comp) <- "age"
GEX.comp
## An object of class Seurat 
## 17900 features across 15536 samples within 1 assay 
## Active assay: RNA (17900 features, 1437 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.comp$age[1:6]
## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1 
##            P30.LPS            P30.LPS            P30.PBS            P30.LPS 
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1 
##            P30.PBS            P30.LPS 
## Levels: P06 P30.PBS P30.LPS
list.b <- list(P30.CTR = c("P30.PBS.CTR","P30.LPS.CTR"),
               P30.MIG = c("P30.PBS.MIG","P30.LPS.MIG"),
               P30.TDT = c("P30.PBS.TDT","P30.LPS.TDT"))
list.b
## $P30.CTR
## [1] "P30.PBS.CTR" "P30.LPS.CTR"
## 
## $P30.MIG
## [1] "P30.PBS.MIG" "P30.LPS.MIG"
## 
## $P30.TDT
## [1] "P30.PBS.TDT" "P30.LPS.TDT"
names(list.b)
## [1] "P30.CTR" "P30.MIG" "P30.TDT"
DEGs.b <- list()

for(nn in names(list.b)){
  
  DEGs.b[[nn]] <- FindAllMarkers(subset(GEX.comp, subset = age %in% c("P30.PBS","P30.LPS") & cnt %in% list.b[[nn]]),
                                 assay = "RNA",
                                 only.pos = T,
                                 min.pct = 0.05,
                                 logfc.threshold = 0.1,
                                 test.use = "MAST")
  
  
}
## Calculating cluster P30.PBS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.PBS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.PBS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
## Calculating cluster P30.LPS
## 
## Done!
## Combining coefficients and standard errors
## Calculating log-fold changes
## Calculating likelihood ratio tests
## Refitting on reduced model...
## 
## Done!
#DEGs.b
names(DEGs.b)
## [1] "P30.CTR" "P30.MIG" "P30.TDT"
# save DEGs' table
df.DEGs.b <- data.frame()

for(nn in names(DEGs.b)){
  df.DEGs.b <- rbind(df.DEGs.b, data.frame(DEGs.b[[nn]], condition = nn))
}

#write.csv(df.DEGs.b, "DEGs.b.subset2.csv")

DEG stat

cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
##   condition cluster up.DEGs
## 1   P30.CTR P30.PBS     425
## 2   P30.CTR P30.LPS     512
## 3   P30.MIG P30.PBS     321
## 4   P30.MIG P30.LPS     391
## 5   P30.TDT P30.PBS     388
## 6   P30.TDT P30.LPS     423
pp.stat.DEG <- list()
df.DEGs.b$cluster <- factor(as.character(df.DEGs.b$cluster),
                            levels = c("P30.PBS","P30.LPS"))
pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]

cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
##   condition cluster up.DEGs
## 1   P30.CTR P30.PBS     425
## 2   P30.CTR P30.LPS     512
## 3   P30.MIG P30.PBS     321
## 4   P30.MIG P30.LPS     391
## 5   P30.TDT P30.PBS     388
## 6   P30.TDT P30.LPS     423
pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]

DEG plot

pp.DEGs.b <- lapply(DEGs.b, function(x){NA})
names(DEGs.b)
## [1] "P30.CTR" "P30.MIG" "P30.TDT"
DEGs.b.combine <- DEGs.b
DEGs.b.combine <- lapply(DEGs.b.combine, function(x){
  x[x$cluster=="P30.PBS","avg_log2FC"] <- -x[x$cluster=="P30.PBS","avg_log2FC"]
  x
})
CTR
pp.DEGs.b$P30.CTR <- DEGs.b.combine$P30.CTR %>%
  mutate(label=ifelse(((p_val_adj < 1e-100 & avg_log2FC<0) | (p_val_adj < 1e-100 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
         padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("P30.PBS"="Skyblue",
                               "P30.LPS"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="CTR, P30.PBS vs P30.LPS") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.CTR

MIG
pp.DEGs.b$P30.MIG <- DEGs.b.combine$P30.MIG %>%
  mutate(label=ifelse(((p_val_adj < 1e-80 & avg_log2FC<0) | (p_val_adj < 1e-80 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
         padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("P30.PBS"="Skyblue",
                               "P30.LPS"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="MIG, P30.PBS vs P30.LPS") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.MIG

TDT
pp.DEGs.b$P30.TDT <- DEGs.b.combine$P30.TDT %>%
  mutate(label=ifelse(((p_val_adj < 1e-80 & avg_log2FC<0) | (p_val_adj < 1e-80 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"P30.PBS",ifelse(p_val_adj<0.05 & avg_log2FC>0,"P30.LPS","None")),
         padj=ifelse(p_val_adj<1e-300,p_val_adj+1e-300,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("P30.PBS"="Skyblue",
                               "P30.LPS"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="TDT, P30.PBS vs P30.LPS") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1,2))
pp.DEGs.b$P30.TDT

signature score

#### 10x data, calculate signature score       
#
## The code below is from Adam Hamber
## 2D scoring by Itay
get_controls <- function(counts, gene.list, verbose=F, control.genes.per.gene=10)
{
    # Itay: "Such scores are inevitably correlated with cell complexity so to avoid 
    # that I subtract a "control" score which is generated by averaging over a control 
    # gene set. Control gene sets are chosen to contain 100 times more genes than the 
    # real gene set (analogous to averaging over 100 control sets of similar size) and 
    # to have the same distribution of population/bulk - based expression levels as the 
    # real gene set, such that they are expected to have the same number of "zeros" and 
    # to eliminate the correlation with complexity."
    # ---------------------------------------------------------------------------------
    # Going to find control points by finding the closest genes in terms of expression level and % of the time we observe it
    if(verbose){cat(sprintf("Finding %s background genes based on similarity to given gene set [%s genes] \n", 
                            control.genes.per.gene*length(gene.list), length(gene.list)))}
    cat("Summarizing data \n")
    summary = data.frame(gene=row.names(counts), mean.expr = Matrix::rowMeans(counts), fract.zero = Matrix::rowMeans(counts==0), stringsAsFactors = F)
    #summary = data.frame(gene=row.names(counts), mean.expr = apply(counts,1,mean), fract.zero = apply(counts==0,1,mean), stringsAsFactors = F)
    summary$mean.expr.s = scale(summary$mean.expr)
    summary$fract.zero.s = scale(summary$fract.zero)
    actual.genes = summary[summary$gene %in% gene.list,]
    background.genes = summary[!summary$gene %in% gene.list,]
    
    #find the 10 closest genes to each cell cycle marker gene and add them to the lists of control genes
    get_closest_genes <- function(i)
    {
        background.genes$dist = sqrt((background.genes$mean.expr.s - actual.genes$mean.expr.s[i])^2 + 
                                         (background.genes$fract.zero.s - actual.genes$fract.zero.s[i])^2)
        ordered = background.genes$gene[order(background.genes$dist)]
        ordered = ordered[!ordered %in% controls] # don't take genes that already appear in the list 
        closest = head(ordered, n=control.genes.per.gene)
        return(closest)
    }
    controls = c();
    
    for (i in 1:length(gene.list)){
        #info(sprintf("Finding %s control genes for %s", control.genes.per.gene, gene.list[i]))
        closest = get_closest_genes(i)
        #info(sprintf("Found %s: ", length(closest)))
        controls = unique(c(controls, closest))
    }
    
    if(verbose){cat(sprintf("Control gene selection complete. %s genes found. \n", length(controls)))}
    #print(controls)
    return(controls)
}

## Define calculate function
calculate_signature_score <- function(count_matrix, gene_list){
    control_gene <- get_controls(counts = count_matrix,
                                 gene.list = gene_list)
    signature_score <- colMeans(count_matrix[gene_list, ], na.rm = TRUE) - 
        colMeans(count_matrix[control_gene, ], na.rm = TRUE)
    return(signature_score)
}

add_geneset_score <- function(obj, geneset, setname){
  score <- calculate_signature_score(as.data.frame(obj@assays[['RNA']]@data),
                                     geneset)
  obj <- AddMetaData(obj,
                     score,
                     setname)
  return(obj)
}
## proc_DEG()
#       to process edgeR result for DEGs' comparison
#              mat_cut is a matrix after advanced filtering, 'like TPM > n in at least one condition'
# 

proc_DEG <- function(deg, p.cut=0.05, FC.cut = 2, padj=TRUE, abs=TRUE, mat_cut=NULL, gene_cut=NULL){
    rownames(deg) <- deg$gene
    
    if(padj==TRUE){
        deg <- deg %>% filter(padj < p.cut)
    }else{
        deg <- deg %>% filter(pvalue < p.cut)
    }
    
    if(abs==TRUE){
        deg <- deg %>% filter(abs(FC) > FC.cut)
    }else if(FC.cut >0){
        deg <- deg %>% filter(FC > FC.cut)
    }else{
        deg <- deg %>% filter(FC < FC.cut)
    }
    
    if(!is.null(mat_cut)){
        deg <- deg[rownames(deg) %in% rownames(mat_cut),]
    }
    if(!is.null(gene_cut)){
        deg <- deg[rownames(deg) %in% gene_cut,]
    }
    return(deg)
}
#saveRDS(GEX.seur, "Spp1tdt.subset2_P30LPS_P30PBS.sort1007.rds")
DEG.list <- list(LPS.TDT= (df.DEGs.a %>% filter(p_val_adj < 0.05 & 
                                               abs(avg_log2FC) > 0.1 & 
                                               pct.1 > 0.05 & contrast == "P30.LPS.TDTvsCTR" &
                                               cluster == "P30.LPS.TDT"))$gene,
                 LPS.CTR= (df.DEGs.a %>% filter(p_val_adj < 0.05 & 
                                               abs(avg_log2FC) > 0.1 & 
                                               pct.1 > 0.05 & contrast == "P30.LPS.TDTvsCTR" &
                                               cluster == "P30.LPS.CTR"))$gene)
DEG.list
## $LPS.TDT
##   [1] "Ly86"          "Xist"          "Cd52"          "Fn1"          
##   [5] "Apoe"          "tdTomato-mid"  "H2-D1"         "H2-K1"        
##   [9] "Cd72"          "Gatm"          "Fth1"          "Ctss"         
##  [13] "Gm21860"       "Ms4a6d"        "B2m"           "Lyz2"         
##  [17] "H2-Q7"         "E230029C05Rik" "Cybb"          "Cd74"         
##  [21] "Ms4a6c"        "Trps1"         "Cp"            "Spp1"         
##  [25] "Ifi204"        "C4b"           "Rps12"         "Fcgr4"        
##  [29] "Rps29"         "Ccl5"          "Ifitm3"        "Ifi211"       
##  [33] "Mpeg1"         "Cd180"         "Saa3"          "Tspo"         
##  [37] "Milr1"         "Ifi207"        "Sash1"         "Chst1"        
##  [41] "Ch25h"         "mt-Nd1"        "H2-Q6"         "Ifi27l2a"     
##  [45] "Axl"           "Npc2"          "Clec4a1"       "AW112010"     
##  [49] "Ctsz"          "Lgals3bp"      "Ccl2"          "Nceh1"        
##  [53] "Rpl39"         "Rplp1"         "Ifi213"        "mt-Atp6"      
##  [57] "Ms4a4c"        "Xdh"           "Malat1"        "Itm2b"        
##  [61] "Lilrb4a"       "Lst1"          "Slfn5"         "Aoah"         
##  [65] "Rps20"         "Bst2"          "Ifitm2"        "Gpr65"        
##  [69] "mt-Co2"        "Zbp1"          "P2rx4"         "Fyb"          
##  [73] "Epb41l3"       "Ptger4"        "Eif2s3y"       "Clec2d"       
##  [77] "H2-Ab1"        "Fcer1g"        "Rps28"         "Lair1"        
##  [81] "Rps27"         "Clec4a3"       "Cpd"           "mt-Cytb"      
##  [85] "Rps24"         "Ifi209"        "Slc6a6"        "Rps3"         
##  [89] "Rps21"         "Hcar2"         "Tor3a"         "Cd48"         
##  [93] "Itgam"         "Iigp1"         "Tmcc3"         "Rpl37a"       
##  [97] "Tlr2"          "Ly9"           "Ifi206"        "Zeb2"         
## [101] "Srgap2"        "AU020206"      "mt-Co3"        "H2-M3"        
## [105] "Ifit3"         "Ptprc"         "Unc93b1"       "Fabp5"        
## [109] "Rsad2"         "Cd14"          "Ccl9"          "Prdx5"        
## [113] "Gm4951"        "Psap"          "Ctsc"          "Mt1"          
## [117] "Tnfsf13b"      "Sp100"         "Rpl23"         "Rasa4"        
## [121] "Anxa5"         "Ccl3"          "Tmem86a"       "Rpl10a"       
## [125] "Ccl6"          "Plgrkt"        "Slc31a2"       "Mbnl1"        
## [129] "Fgl2"          "Oasl2"        
## 
## $LPS.CTR
##   [1] "tdTomato-head" "Hspa8"         "Dnaja1"        "Hsp90ab1"     
##   [5] "Hsp90aa1"      "Cd164"         "Fcrls"         "Tmem119"      
##   [9] "Hspa1a"        "Hspa1b"        "Hspa5"         "Hsph1"        
##  [13] "Cacybp"        "Hspd1"         "Hspe1"         "G530011O06Rik"
##  [17] "Pmepa1"        "Fkbp4"         "Lrba"          "Sparc"        
##  [21] "Tuba1a"        "Basp1"         "Stip1"         "Mat2a"        
##  [25] "Gt(ROSA)26Sor" "Abi3"          "Chordc1"       "Slc2a5"       
##  [29] "Manf"          "Ubb"           "Cct3"          "Cebpd"        
##  [33] "Eif1"          "Cd34"          "P4ha1"         "Ptges3"       
##  [37] "Ubc"           "Cct7"          "Cct6a"         "Eif5"         
##  [41] "Dnajb1"        "Srm"           "Fam102b"       "Eif2s2"       
##  [45] "P2ry12"        "Cebpb"         "Calm2"         "Set"          
##  [49] "Ogfrl1"        "C5ar1"         "Pten"          "Phgdh"        
##  [53] "Bin1"          "Eif4a1"        "Cct5"          "Ppp1r14b"     
##  [57] "Rab3il1"       "Ahsa2"         "Rps6ka1"       "Jam3"         
##  [61] "Denr"          "St13"          "Ifi27"         "Serbp1"       
##  [65] "Irf2bp2"       "Tcf7l2"        "Mgat2"         "Snx18"        
##  [69] "Npl"           "Top1"          "Nars"          "Fchsd2"       
##  [73] "Entpd1"        "Ptma"          "Chsy1"         "Ckb"          
##  [77] "Mrfap1"        "Arpc2"         "Picalm"        "Cct4"         
##  [81] "Fscn1"         "Bag3"          "Mcrip1"        "Fxr1"         
##  [85] "Selenos"       "Orai3"         "Txndc11"       "Ywhae"        
##  [89] "Cdk6"          "Ltc4s"         "Spsb1"         "Pcbp1"        
##  [93] "Creld2"        "Mertk"         "Gbp7"          "Ranbp1"       
##  [97] "Nifk"          "Scoc"          "Klf3"          "Magoh"        
## [101] "Ran"           "Hmgb1"         "Npm1"          "Tspan7"       
## [105] "Bank1"         "Tubb5"         "Nrip1"         "Arhgdia"      
## [109] "Plxdc2"        "Tmem65"        "Spcs2"         "Cct2"         
## [113] "Pnp"           "Ptgs1"         "Tmem204"       "Pafah1b1"     
## [117] "Lsm2"          "Ryk"           "Cib1"          "Gas6"         
## [121] "Slc39a10"      "Washc4"        "Bin2"          "Fli1"         
## [125] "Sel1l"         "Tcea1"         "Adgrg1"        "Csk"          
## [129] "Idh2"          "Crebbp"        "Ipo7"          "Tuba1b"       
## [133] "Usp10"         "Morf4l2"       "Atp2a2"        "Clta"         
## [137] "Spint1"        "Mrps6"         "Gcnt1"         "G3bp1"        
## [141] "Il1a"          "Sdf2l1"        "Mrps12"        "Rab14"        
## [145] "Nfam1"         "Dynll1"        "Vapa"          "Cd300a"       
## [149] "4833439L19Rik" "Ppa1"          "Eif1b"         "Ifrd1"        
## [153] "B3galt5"       "Pa2g4"         "Rps19bp1"      "Rapgef5"      
## [157] "Klhl9"         "Eif3c"         "Ywhah"         "Erh"          
## [161] "Pip4k2a"       "Pak2"          "Sox4"          "Ebna1bp2"     
## [165] "Cnbp"          "Tjp1"          "Tmem173"       "Alox5ap"      
## [169] "Ppp2ca"        "Tagln2"        "Rbpj"          "Mef2a"        
## [173] "Rpl36al"       "Mak16"         "Rangap1"       "Eif3j1"       
## [177] "Mtus1"         "Fam167b"       "Rbm8a"         "Slc25a5"      
## [181] "Eif4e"         "Ncl"           "Scaf11"        "Nudc"         
## [185] "Ahsa1"         "Hnrnpu"        "Smndc1"        "Snx6"         
## [189] "Fez2"          "Mex3c"         "Sall1"         "Rrp1b"        
## [193] "Zfr"           "Eif4g2"        "Cct8"          "Yy1"          
## [197] "Gtf2i"         "Tmed7"         "Mylip"         "Cltc"         
## [201] "Rexo2"         "Comt"          "Gmps"          "Dctpp1"       
## [205] "Gngt2"         "C1qbp"         "Ubxn4"         "Nolc1"        
## [209] "Khk"
lapply(DEG.list, length)
## $LPS.TDT
## [1] 130
## 
## $LPS.CTR
## [1] 209
GEX.seur <- add_geneset_score(GEX.seur, DEG.list$LPS.TDT, "LPS.TDT.up")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DEG.list$LPS.CTR, "LPS.TDT.dn")
## Summarizing data
ppnew.2.v1 <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn"), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
ppnew.2.v1

ppnew.2.v1a <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn")[1], 
                      same.y.lims = T,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT"),
                                                  c("P30.PBS.TDT","P30.LPS.TDT")),
                               label.y = c(0.55,0.75,0.95,0.55,0.75,0.95,1.1),
                               size=3
                               )
ppnew.2.v1a

ppnew.2.v1b <- VlnPlot(GEX.seur, features = c("LPS.TDT.up","LPS.TDT.dn")[2], 
                      same.y.lims = T,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & ylim(c(-0.23,0.57)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT")#,
                                                  #c("P30.PBS.TDT","P30.LPS.TDT")
                                                  ),
                               label.y = c(0.25,0.3,0.4,0.35,0.4,0.5,0.55),
                               size=3
                               )
ppnew.2.v1b

DotPlot(GEX.seur, features = c("LPS.TDT.dn","LPS.TDT.up"), group.by = "cnt", cols = c("midnightblue","darkorange1"), scale = F, scale.by = "size", scale.min = 1, scale.max = 1, dot.scale = 12)  + 
  coord_flip() + labs(title="mean signature score") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) 

GEX.seur@meta.data[,c("cnt","LPS.TDT.dn","LPS.TDT.up")] %>% 
  reshape2::melt("cnt") %>%
  ggplot(aes(x = cnt, y = variable, color=value)) +
  geom_point(size=12) +
  scale_color_gradientn(colors = c("midnightblue","darkorange1")) +
  theme_classic()  + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

new clusters

GEX.seur$new_clusters <- as.character(GEX.seur$seurat_clusters)

GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(1)] <- "C1"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(2,8,7,14)] <- "C2"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(9)] <- "C3"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(12)] <- "C4"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(3)] <- "C5"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(5)] <- "C6"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(0,4)] <- "C7"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(6)] <- "C8"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(11)] <- "C9"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(13)] <- "C10"
GEX.seur$new_clusters[GEX.seur$new_clusters %in% c(10)] <- "C11"

GEX.seur$new_clusters <- factor(GEX.seur$new_clusters,
                                levels = paste0("C",1:11))
head(GEX.seur@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGGTCGACA-1    Spp1tdt       5001         2062   2.699460  12.437512
## AAACCCAAGGTTCTTG-1    Spp1tdt       6501         2457   2.553453  12.305799
## AAACCCACAAACTCGT-1    Spp1tdt       4188         1841   2.960840   8.572111
## AAACCCACACGCGGTT-1    Spp1tdt       4524         1913   1.215738   9.858532
## AAACCCACATCCGTGG-1    Spp1tdt       2337         1225   3.979461   7.659392
## AAACCCAGTATCGAAA-1    Spp1tdt       4318         1955   2.292728  11.741547
##                         FB.info       S.Score   G2M.Score Phase RNA_snn_res.1.5
## AAACCCAAGGTCGACA-1  P30.LPS.CTR  0.0061184939 -0.06523366     S               0
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT1  0.0002076412 -0.07102627     S               0
## AAACCCACAAACTCGT-1  P30.PBS.MIG  0.0192137320  0.06031440   G2M               7
## AAACCCACACGCGGTT-1 P30.LPS.TDT2 -0.0535299003 -0.05464656    G1               0
## AAACCCACATCCGTGG-1  P30.PBS.TDT  0.0489756368 -0.04211671     S               2
## AAACCCAGTATCGAAA-1  P30.LPS.MIG -0.0314645626  0.02183367   G2M               0
##                    seurat_clusters pANN_0.25_0.005_1191 DoubletFinder0.05
## AAACCCAAGGTCGACA-1               0          0.000000000           Singlet
## AAACCCAAGGTTCTTG-1               0          0.044025157           Singlet
## AAACCCACAAACTCGT-1               7          0.006289308           Singlet
## AAACCCACACGCGGTT-1               0          0.012578616           Singlet
## AAACCCACATCCGTGG-1               2          0.000000000           Singlet
## AAACCCAGTATCGAAA-1               0          0.006289308           Singlet
##                    pANN_0.25_0.005_2383 DoubletFinder0.1 sort_clusters     age
## AAACCCAAGGTCGACA-1          0.000000000          Singlet             0 P30.LPS
## AAACCCAAGGTTCTTG-1          0.044025157          Singlet             0 P30.LPS
## AAACCCACAAACTCGT-1          0.006289308          Singlet             7 P30.PBS
## AAACCCACACGCGGTT-1          0.012578616          Singlet             0 P30.LPS
## AAACCCACATCCGTGG-1          0.000000000          Singlet             2 P30.PBS
## AAACCCAGTATCGAAA-1          0.006289308          Singlet             0 P30.LPS
##                            cnt   preAnno RNA_snn_res.1 RNA_snn_res.0.8
## AAACCCAAGGTCGACA-1 P30.LPS.CTR Microglia             2               2
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT Microglia             3               4
## AAACCCACAAACTCGT-1 P30.PBS.MIG Microglia             7               0
## AAACCCACACGCGGTT-1 P30.LPS.TDT Microglia             2               2
## AAACCCACATCCGTGG-1 P30.PBS.TDT Microglia             0               0
## AAACCCAGTATCGAAA-1 P30.LPS.MIG Microglia             2               2
##                    RNA_snn_res.0.6       FB.new RNA_snn_res.1.2 RNA_snn_res.1.3
## AAACCCAAGGTCGACA-1               0  P30.LPS.CTR               2               2
## AAACCCAAGGTTCTTG-1               0 P30.LPS.TDT1               3               4
## AAACCCACAAACTCGT-1               1  P30.PBS.MIG               7               7
## AAACCCACACGCGGTT-1               0 P30.LPS.TDT2               3               4
## AAACCCACATCCGTGG-1               1  P30.PBS.TDT               1               1
## AAACCCAGTATCGAAA-1               0  P30.LPS.MIG               2               2
##                    RNA_snn_res.1.4 LPS.TDT.up  LPS.TDT.dn new_clusters    age1
## AAACCCAAGGTCGACA-1               2  0.2671683  0.24670044           C7 P30.LPS
## AAACCCAAGGTTCTTG-1               2  0.2190005  0.21501295           C7 P30.LPS
## AAACCCACAAACTCGT-1               7  0.3136800  0.02545367           C2 P30.PBS
## AAACCCACACGCGGTT-1               2  0.2974201  0.09338095           C7 P30.LPS
## AAACCCACATCCGTGG-1               1  0.2758054 -0.03406151           C2 P30.PBS
## AAACCCAGTATCGAAA-1               2  0.2300058  0.11176488           C7 P30.LPS
##                       cluster_cnt DAM.sig_top50 DAM.sig_top100 DAM.sig_top250
## AAACCCAAGGTCGACA-1 C7_P30.LPS.CTR    0.12647123     0.09610184     0.09546331
## AAACCCAAGGTTCTTG-1 C7_P30.LPS.TDT    0.04136910     0.02314892     0.06780191
## AAACCCACAAACTCGT-1 C2_P30.PBS.MIG   -0.03826143    -0.02304032    -0.01821938
## AAACCCACACGCGGTT-1 C7_P30.LPS.TDT    0.17479143     0.26609309     0.19238189
## AAACCCACATCCGTGG-1 C2_P30.PBS.TDT    0.06229684    -0.06607385    -0.05013960
## AAACCCAGTATCGAAA-1 C7_P30.LPS.MIG    0.19221428     0.07223828     0.07965503
##                    DAM.sig_top500 LPS.sig_score
## AAACCCAAGGTCGACA-1      0.2497200     0.5988063
## AAACCCAAGGTTCTTG-1      0.2634468     0.6298800
## AAACCCACAAACTCGT-1      0.1498645     0.2308747
## AAACCCACACGCGGTT-1      0.3066293     0.6339013
## AAACCCACATCCGTGG-1      0.1023302     0.1488521
## AAACCCAGTATCGAAA-1      0.2186862     0.4894010
#
GEX.seur$age1 <- factor(as.character(GEX.seur$age),
                        levels = c("P30.PBS","P30.LPS"))

#
GEX.seur$FB.new <- factor(as.character(GEX.seur$FB.info),
                          levels = rev(levels(GEX.seur$FB.info))[c(8:10,4,5,7,6)])

color.FBnew <- color.FB2[c(3:1,7,6,4,5)]
head(GEX.seur$cnt)
## AAACCCAAGGTCGACA-1 AAACCCAAGGTTCTTG-1 AAACCCACAAACTCGT-1 AAACCCACACGCGGTT-1 
##        P30.LPS.CTR        P30.LPS.TDT        P30.PBS.MIG        P30.LPS.TDT 
## AAACCCACATCCGTGG-1 AAACCCAGTATCGAAA-1 
##        P30.PBS.TDT        P30.LPS.MIG 
## 6 Levels: P30.PBS.CTR P30.PBS.MIG P30.PBS.TDT P30.LPS.CTR ... P30.LPS.TDT
DotPlot(GEX.seur, features = c(markers.mig,"Ifit1","Ifit2"), group.by = "new_clusters",
                         cols = c("midnightblue","darkorange1")) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1))#+ scale_y_discrete(limits=rev)

new markers

# find markers for every cluster compared to all remaining cells, report only the positive ones
Idents(GEX.seur) <- "new_clusters"

#GEX.markers.new <- FindAllMarkers(GEX.seur, only.pos = TRUE, 
#                                  min.pct = 0.05,
#                                  test.use = "MAST",
#                                  logfc.threshold = 0.25)
GEX.markers.new <- read.table("sc10x_LYN.0921.marker.subset2_new.csv", header = TRUE, sep = ",")
GEX.markers.new %>% group_by(cluster) %>% top_n(n = 8, wt = avg_log2FC)
## # A tibble: 88 x 7
## # Groups:   cluster [11]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene   
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <chr>   <chr>  
##  1 0              1.05  0.999 0.973 0         C1      P2ry12 
##  2 0              0.965 0.998 0.921 0         C1      Gpr34  
##  3 0              0.812 0.979 0.897 0         C1      Tgfbr1 
##  4 2.09e-269      0.751 0.966 0.857 3.73e-265 C1      Arhgap5
##  5 1.65e-265      0.731 0.973 0.903 2.96e-261 C1      Maf    
##  6 4.40e-250      0.748 0.925 0.764 7.88e-246 C1      Rhob   
##  7 9.96e-142      0.717 0.511 0.276 1.78e-137 C1      Slc40a1
##  8 3.05e-121      0.717 0.566 0.329 5.46e-117 C1      Sox4   
##  9 0              1.03  0.993 0.91  0         C2      Gpr34  
## 10 0              1.01  1     0.968 0         C2      P2ry12 
## # ... with 78 more rows
markers.new <- GEX.markers.new
markers.new$cluster <- factor(as.character(GEX.markers.new$cluster),
                               levels = levels(GEX.seur$new_clusters))

markers.new_t60 <- (markers.new %>% group_by(cluster) %>% 
                  filter(pct.1>0.1 & gene %in% grep("Rps|Rpl|Gm|Rik$|Xist|Tsix",markers.new$gene,invert = T,value = T)) %>%
                   top_n(n = 60, wt = avg_log2FC) %>%
                   ungroup() %>%
  #arrange(desc(avg_log2FC*pct.1),gene) %>%
    arrange(desc(pct.1),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[1:64]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[65:128]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[129:192]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[193:256]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[257:320]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[321:384]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

DotPlot(GEX.seur, group.by = "new_clusters", features = rev(markers.new_t60[385:436]))  + coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6))

plot1

umap in clusters

pp.umap.1a <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "seurat_clusters", raster = F, pt.size = 0.3) +
           DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "FB.new", cols = color.FBnew, raster = F, pt.size = 0.3)
pp.umap.1a

pp.umap.1b <- DimPlot(GEX.seur, reduction = "umap", label = T, group.by = "new_clusters", repel = F, raster = F, pt.size = 0.3) +
           DimPlot(GEX.seur, reduction = "umap", label = F, group.by = "cnt", cols = color.cnt, raster = F, pt.size = 0.3)
pp.umap.1b

plot2

pp.umap.2a <- DimPlot(GEX.seur, label = F,  repel = F, reduction = 'umap', group.by = "cnt", split.by = "cnt", raster = F, pt.size = 0.3,
        ncol = 3, cols = color.cnt)
pp.umap.2a

pp.umap.2b <- DimPlot(GEX.seur, label = F,repel = F, reduction = 'umap', group.by = "cnt", split.by = "age", 
                     raster = F, pt.size = 0.3,
                     ncol = 3, cols = color.cnt)
pp.umap.2b

pp.umap.2c <- DimPlot(GEX.seur, label = F,  repel = F, reduction = 'umap', group.by = "new_clusters", split.by = "cnt", raster = F, pt.size = 0.3,
        ncol = 3)
pp.umap.2c

pp.umap.2d <- DimPlot(GEX.seur, label = F,repel = F, reduction = 'umap', group.by = "new_clusters", split.by = "age", 
                     raster = F, pt.size = 0.3,
                     ncol = 3)
pp.umap.2d

plot3

composition

pp.comp.3a <- cowplot::plot_grid(
pheatmap::pheatmap(table(FB.info=GEX.seur$cnt,
      clusters=GEX.seur$new_clusters),
                   main = "Cell Count",
                   gaps_row = c(3),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.0f", legend = F, silent = T)$gtable,

pheatmap::pheatmap(100*rowRatio(table(FB.info=GEX.seur$cnt,
                                clusters=GEX.seur$new_clusters)),
                   main = "Cell Ratio",
                   gaps_row = c(3),
      #gaps_col = c(12,13,14),
                   cluster_rows = F, cluster_cols = F, display_numbers = T, number_format = "%.1f", legend = F, silent =T)$gtable,
ncol = 1)
pp.comp.3a

stat_sort <- GEX.seur@meta.data[,c("cnt","new_clusters")]
stat_sort.s <- stat_sort %>%
  group_by(cnt,new_clusters) %>%
  summarise(count=n()) %>%
  mutate(prop= count/sum(count)) %>%
  ungroup()
## `summarise()` has grouped output by 'cnt'. You can override using the `.groups`
## argument.
dim(stat_sort.s)
## [1] 55  4
#stat_sort.s
pp.comp.3b <- stat_sort.s %>% 
  ggplot(aes(x = new_clusters, y = 100*prop, color=cnt)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  theme_classic(base_size = 15) + 
  scale_color_manual(values = color.FBnew, name="") +
  labs(title="Cell Composition", y="Proportion") + 
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6)) +
  
  geom_point(aes(x = new_clusters, y = 100*prop, color=cnt),
             position = position_dodge(0.75),
             shape=16,alpha=0.75,size=1.75,
             stroke=0.15, show.legend = F)
pp.comp.3b

plot4

features4 <- c("Saa3","Il1a","Ccl3","Ccl4",
               "Ccl5","Il1b","Tnf","Ccl2",
               "Cxcl10")
pp.feat.4 <- FeaturePlot(GEX.seur, 
                         features = features4, raster = T, pt.size = 3.5,
                         ncol = 4)
pp.feat.4

plot5

final.markers <- c("Apoe","Ms4a6d","Ms4a6c","Ctss",
                   "Ctsc","Ctsd","Ch25h","Ccl2",
                   "Ccl3","Ccl4","Ccl5","Cxcl10",
                   "Il1a","Il1b","Tnf","Saa3")
final.markers
##  [1] "Apoe"   "Ms4a6d" "Ms4a6c" "Ctss"   "Ctsc"   "Ctsd"   "Ch25h"  "Ccl2"  
##  [9] "Ccl3"   "Ccl4"   "Ccl5"   "Cxcl10" "Il1a"   "Il1b"   "Tnf"    "Saa3"
length(final.markers)
## [1] 16
sum(final.markers %in% markers.new$gene)
## [1] 16
final.markers.sort <- (markers.new %>% group_by(cluster) %>%
                         filter(gene %in% final.markers) %>%
                         ungroup() %>%
                         arrange(desc(avg_log2FC*pct.1),gene) %>%
                         distinct(gene, .keep_all = TRUE) %>%
                         arrange(cluster, p_val_adj))$gene
final.markers.sort
##  [1] "Ctsd"   "Il1a"   "Saa3"   "Ms4a6c" "Ch25h"  "Ms4a6d" "Ctsc"   "Ccl5"  
##  [9] "Ccl3"   "Tnf"    "Il1b"   "Cxcl10" "Ccl2"   "Ccl4"   "Apoe"   "Ctss"
final.markers.new <- final.markers.sort
mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
pp.dot.5a1 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.5a1

pp.dot.5a2 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")  +
  scale_color_gradientn(colours  = rev(mapal))
pp.dot.5a2

pp.dot.5b1 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions")
pp.dot.5b1

pp.dot.5b2 <- DotPlot(GEX.seur, features = final.markers.new, group.by = "cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions")  +
  scale_color_gradientn(colours  = rev(mapal))
pp.dot.5b2

levels.cluster_cnt <- as.vector(unlist(sapply(levels(GEX.seur$new_clusters),function(x){
                         paste0(x,"_",levels(GEX.seur$cnt))
                      })))
levels.cluster_cnt
##  [1] "C1_P30.PBS.CTR"  "C1_P30.PBS.MIG"  "C1_P30.PBS.TDT"  "C1_P30.LPS.CTR" 
##  [5] "C1_P30.LPS.MIG"  "C1_P30.LPS.TDT"  "C2_P30.PBS.CTR"  "C2_P30.PBS.MIG" 
##  [9] "C2_P30.PBS.TDT"  "C2_P30.LPS.CTR"  "C2_P30.LPS.MIG"  "C2_P30.LPS.TDT" 
## [13] "C3_P30.PBS.CTR"  "C3_P30.PBS.MIG"  "C3_P30.PBS.TDT"  "C3_P30.LPS.CTR" 
## [17] "C3_P30.LPS.MIG"  "C3_P30.LPS.TDT"  "C4_P30.PBS.CTR"  "C4_P30.PBS.MIG" 
## [21] "C4_P30.PBS.TDT"  "C4_P30.LPS.CTR"  "C4_P30.LPS.MIG"  "C4_P30.LPS.TDT" 
## [25] "C5_P30.PBS.CTR"  "C5_P30.PBS.MIG"  "C5_P30.PBS.TDT"  "C5_P30.LPS.CTR" 
## [29] "C5_P30.LPS.MIG"  "C5_P30.LPS.TDT"  "C6_P30.PBS.CTR"  "C6_P30.PBS.MIG" 
## [33] "C6_P30.PBS.TDT"  "C6_P30.LPS.CTR"  "C6_P30.LPS.MIG"  "C6_P30.LPS.TDT" 
## [37] "C7_P30.PBS.CTR"  "C7_P30.PBS.MIG"  "C7_P30.PBS.TDT"  "C7_P30.LPS.CTR" 
## [41] "C7_P30.LPS.MIG"  "C7_P30.LPS.TDT"  "C8_P30.PBS.CTR"  "C8_P30.PBS.MIG" 
## [45] "C8_P30.PBS.TDT"  "C8_P30.LPS.CTR"  "C8_P30.LPS.MIG"  "C8_P30.LPS.TDT" 
## [49] "C9_P30.PBS.CTR"  "C9_P30.PBS.MIG"  "C9_P30.PBS.TDT"  "C9_P30.LPS.CTR" 
## [53] "C9_P30.LPS.MIG"  "C9_P30.LPS.TDT"  "C10_P30.PBS.CTR" "C10_P30.PBS.MIG"
## [57] "C10_P30.PBS.TDT" "C10_P30.LPS.CTR" "C10_P30.LPS.MIG" "C10_P30.LPS.TDT"
## [61] "C11_P30.PBS.CTR" "C11_P30.PBS.MIG" "C11_P30.PBS.TDT" "C11_P30.LPS.CTR"
## [65] "C11_P30.LPS.MIG" "C11_P30.LPS.TDT"
GEX.seur$cluster_cnt <- factor(paste0(as.character(GEX.seur$new_clusters),
                                      "_",
                                      as.character(GEX.seur$cnt)),
                               levels = levels.cluster_cnt)
pp.dot.5c1a <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                       features = final.markers.new, group.by = "cluster_cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions")
pp.dot.5c1a

pp.dot.5c2a <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                       features = final.markers.new, group.by = "cluster_cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions") +
  scale_color_gradientn(colours  = rev(mapal))
pp.dot.5c2a

pp.dot.5c1b <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                       features = final.markers.new, group.by = "cluster_cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions")
pp.dot.5c1b

pp.dot.5c2b <- DotPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                       features = final.markers.new, group.by = "cluster_cnt"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters_Conditions") +
  scale_color_gradientn(colours  = rev(mapal))
pp.dot.5c2b

plot6

violin using same markers in plot5

pp.violin.6a1 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters", 
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6a1

pp.violin.6a2 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "new_clusters", 
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6a2

pp.violin.6b1 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6b1

pp.violin.6b2 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6b2

pp.violin.6b3 <- VlnPlot(GEX.seur, features = final.markers.new, group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,7.2)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.TDT","P30.PBS.TDT")),
                               label.y = c(6.3,6.8,7.3,6.3,6.8,7.3,7.8)-1,
                               size=1.75
                               )
pp.violin.6b3

pp.violin.6c1a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6c1a

pp.violin.6c1b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.violin.6c1b

pp.violin.6c2a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6c2a

pp.violin.6c2b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.6c2b

contrast.a <- list(c("P30.PBS.CTR","P30.PBS.MIG"),
                   c("P30.PBS.MIG","P30.PBS.TDT"),
                   c("P30.PBS.CTR","P30.PBS.TDT"))

plot.contrast.a <- list()

for(XXX in paste0("C",1:3)){
  
  
  plot.contrast.a <- c(plot.contrast.a,
                       lapply(contrast.a, function(x){
                         paste0(XXX,
                                "_",
                                x)
                       }))
}

plot.labely.a <- rep(c(6.3,6.8,7.3),3)
#
length(plot.contrast.a)
## [1] 9
plot.contrast.a[1:6]
## [[1]]
## [1] "C1_P30.PBS.CTR" "C1_P30.PBS.MIG"
## 
## [[2]]
## [1] "C1_P30.PBS.MIG" "C1_P30.PBS.TDT"
## 
## [[3]]
## [1] "C1_P30.PBS.CTR" "C1_P30.PBS.TDT"
## 
## [[4]]
## [1] "C2_P30.PBS.CTR" "C2_P30.PBS.MIG"
## 
## [[5]]
## [1] "C2_P30.PBS.MIG" "C2_P30.PBS.TDT"
## 
## [[6]]
## [1] "C2_P30.PBS.CTR" "C2_P30.PBS.TDT"
contrast.b <- list(c("P30.LPS.CTR","P30.LPS.MIG"),
                   c("P30.LPS.MIG","P30.LPS.TDT"),
                   c("P30.LPS.CTR","P30.LPS.TDT"))

plot.contrast.b <- list()

for(XXX in paste0("C",4:11)){
  
  
  plot.contrast.b <- c(plot.contrast.b,
                       lapply(contrast.b, function(x){
                         paste0(XXX,
                                "_",
                                x)
                       }))
}

plot.labely.b <- rep(c(6.3,6.8,7.3),8)
#
length(plot.contrast.b)
## [1] 24
plot.contrast.b[1:6]
## [[1]]
## [1] "C4_P30.LPS.CTR" "C4_P30.LPS.MIG"
## 
## [[2]]
## [1] "C4_P30.LPS.MIG" "C4_P30.LPS.TDT"
## 
## [[3]]
## [1] "C4_P30.LPS.CTR" "C4_P30.LPS.TDT"
## 
## [[4]]
## [1] "C5_P30.LPS.CTR" "C5_P30.LPS.MIG"
## 
## [[5]]
## [1] "C5_P30.LPS.MIG" "C5_P30.LPS.TDT"
## 
## [[6]]
## [1] "C5_P30.LPS.CTR" "C5_P30.LPS.TDT"
pp.violin.6c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.FBnew[1:3],4),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)& ylim(c(0,6.8)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.a,
                               label.y = plot.labely.a-1,
                               size=1.5
                               ) 
pp.violin.6c3a

pp.violin.6c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                          features = final.markers.new, group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)& ylim(c(0,6.8)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.b,
                               label.y = plot.labely.b-1,
                               size=1.75
                               )
pp.violin.6c3b

plot7

count DEG number, using different cutoffs

cut.padj = 0.05
cut.log2FC = 0.1   
  
cut.pct1 = 0.05

stat.a1 <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat.a1
##           contrast     cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR     209
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT     130
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG     187
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT     164
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR      28
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT      35
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG      21
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT       5
#df.DEGs.a
pp.stat.DEG[["a1"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["a1"]]

cut.padj = 0.05
cut.log2FC = 0.2   
  
cut.pct1 = 0.05

stat.a2 <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat.a2
##           contrast     cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR      30
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT      53
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG      44
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT      94
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR       9
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT      17
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG       9
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT       4
pp.stat.DEG[["a2"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["a2"]]

cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

stat.a3 <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
stat.a3
##           contrast     cluster up.DEGs
## 1 P30.LPS.TDTvsCTR P30.LPS.CTR       8
## 2 P30.LPS.TDTvsCTR P30.LPS.TDT      29
## 3 P30.LPS.TDTvsMIG P30.LPS.MIG      17
## 4 P30.LPS.TDTvsMIG P30.LPS.TDT      51
## 5 P30.PBS.TDTvsCTR P30.PBS.CTR       2
## 6 P30.PBS.TDTvsCTR P30.PBS.TDT       7
## 7 P30.PBS.TDTvsMIG P30.PBS.MIG       3
## 8 P30.PBS.TDTvsMIG P30.PBS.TDT       1
pp.stat.DEG[["a3"]] <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=contrast, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = color.cnt[c(4:6,1:3)], name="") +
  labs(title=paste0("up.DEGs stat, pct.1>0.05, padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["a3"]]

DEG list

cut.padj = 0.05
cut.log2FC = 0.1   
  
cut.pct1 = 0.05

stat.a1.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  #summarise(up.DEGs = n()) %>% 
  as.data.frame()
dim(stat.a1.df)
## [1] 779   8
stat.a1.df[1:6,]
##          p_val avg_log2FC pct.1 pct.2    p_val_adj     cluster          gene
## 1 7.291059e-54  0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 5.510551e-20  0.1877754 0.998 0.990 9.863886e-16 P30.PBS.CTR         Fcrls
## 3 7.112810e-19  0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR          Lrba
## 4 2.383831e-13  0.1595905 0.999 0.999 4.267058e-09 P30.PBS.CTR        P2ry12
## 5 7.947888e-13  0.1056330 1.000 1.000 1.422672e-08 P30.PBS.CTR          Actb
## 6 1.076088e-12  0.2209830 0.892 0.837 1.926198e-08 P30.PBS.CTR         Cd164
##           contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTR
cut.padj = 0.05
cut.log2FC = 0.2   
  
cut.pct1 = 0.05

stat.a2.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  #summarise(up.DEGs = n()) %>% 
  as.data.frame()
dim(stat.a2.df)
## [1] 260   8
stat.a2.df[1:6,]
##          p_val avg_log2FC pct.1 pct.2    p_val_adj     cluster          gene
## 1 7.291059e-54  0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 7.112810e-19  0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR          Lrba
## 3 1.076088e-12  0.2209830 0.892 0.837 1.926198e-08 P30.PBS.CTR         Cd164
## 4 2.513301e-12  0.2151293 0.892 0.828 4.498809e-08 P30.PBS.CTR        Pmepa1
## 5 3.785970e-12  0.2683232 0.287 0.192 6.776886e-08 P30.PBS.CTR         Bank1
## 6 2.500566e-11  0.2553055 0.137 0.071 4.476013e-07 P30.PBS.CTR         Upk1b
##           contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTR
cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

stat.a3.df <- df.DEGs.a %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(contrast,cluster) %>%
  #summarise(up.DEGs = n()) %>% 
  as.data.frame()
dim(stat.a3.df)
## [1] 118   8
stat.a3.df[1:6,]
##          p_val avg_log2FC pct.1 pct.2    p_val_adj     cluster          gene
## 1 7.291059e-54  0.5042723 0.174 0.033 1.305100e-49 P30.PBS.CTR tdTomato-head
## 2 7.112810e-19  0.3197571 0.447 0.310 1.273193e-14 P30.PBS.CTR          Lrba
## 3 3.177818e-55  0.6194995 0.827 0.722 5.688295e-51 P30.PBS.TDT         H2-D1
## 4 1.418235e-47  0.5638234 0.819 0.743 2.538641e-43 P30.PBS.TDT         H2-K1
## 5 1.868002e-28  0.7477478 0.966 0.957 3.343723e-24 P30.PBS.TDT          Apoe
## 6 1.724448e-20  0.5015367 0.495 0.400 3.086762e-16 P30.PBS.TDT          Lyz2
##           contrast
## 1 P30.PBS.TDTvsCTR
## 2 P30.PBS.TDTvsCTR
## 3 P30.PBS.TDTvsCTR
## 4 P30.PBS.TDTvsCTR
## 5 P30.PBS.TDTvsCTR
## 6 P30.PBS.TDTvsCTR
#
stat.a1.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
                     P30.LPS.TDTvsCTR.TDTup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
                     P30.LPS.TDTvsMIG.MIGup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
                     P30.LPS.TDTvsMIG.TDTup = (stat.a1.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
                     P30.PBS.TDTvsCTR.CTRup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
                     P30.PBS.TDTvsCTR.TDTup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
                     P30.PBS.TDTvsMIG.MIGup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
                     P30.PBS.TDTvsMIG.TDTup = (stat.a1.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)
stat.a2.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
                     P30.LPS.TDTvsCTR.TDTup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
                     P30.LPS.TDTvsMIG.MIGup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
                     P30.LPS.TDTvsMIG.TDTup = (stat.a2.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
                     P30.PBS.TDTvsCTR.CTRup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
                     P30.PBS.TDTvsCTR.TDTup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
                     P30.PBS.TDTvsMIG.MIGup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
                     P30.PBS.TDTvsMIG.TDTup = (stat.a2.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)
stat.a3.list <- list(P30.LPS.TDTvsCTR.CTRup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.CTR"))$gene,
                     P30.LPS.TDTvsCTR.TDTup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsCTR" & cluster == "P30.LPS.TDT"))$gene,
                     P30.LPS.TDTvsMIG.MIGup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.MIG"))$gene,
                     P30.LPS.TDTvsMIG.TDTup = (stat.a3.df %>% filter(contrast == "P30.LPS.TDTvsMIG" & cluster == "P30.LPS.TDT"))$gene,
                     P30.PBS.TDTvsCTR.CTRup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.CTR"))$gene,
                     P30.PBS.TDTvsCTR.TDTup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsCTR" & cluster == "P30.PBS.TDT"))$gene,
                     P30.PBS.TDTvsMIG.MIGup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.MIG"))$gene,
                     P30.PBS.TDTvsMIG.TDTup = (stat.a3.df %>% filter(contrast == "P30.PBS.TDTvsMIG" & cluster == "P30.PBS.TDT"))$gene)
lapply(stat.a1.list,length)
## $P30.LPS.TDTvsCTR.CTRup
## [1] 209
## 
## $P30.LPS.TDTvsCTR.TDTup
## [1] 130
## 
## $P30.LPS.TDTvsMIG.MIGup
## [1] 187
## 
## $P30.LPS.TDTvsMIG.TDTup
## [1] 164
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 28
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 35
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 21
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 5
lapply(stat.a2.list,length)
## $P30.LPS.TDTvsCTR.CTRup
## [1] 30
## 
## $P30.LPS.TDTvsCTR.TDTup
## [1] 53
## 
## $P30.LPS.TDTvsMIG.MIGup
## [1] 44
## 
## $P30.LPS.TDTvsMIG.TDTup
## [1] 94
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 9
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 17
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 9
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 4
lapply(stat.a3.list,length)
## $P30.LPS.TDTvsCTR.CTRup
## [1] 8
## 
## $P30.LPS.TDTvsCTR.TDTup
## [1] 29
## 
## $P30.LPS.TDTvsMIG.MIGup
## [1] 17
## 
## $P30.LPS.TDTvsMIG.TDTup
## [1] 51
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 2
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 7
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 3
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 1
load subset1 DEGs.a list
#
stat.s1.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsCTR.CTRup.txt"))),
                     P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsCTR.TDTup.txt"))),
                     P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsMIG.MIGup.txt"))),
                     P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P06.TDTvsMIG.TDTup.txt"))),
                     P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsCTR.CTRup.txt"))),
                     P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsCTR.TDTup.txt"))),
                     P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsMIG.MIGup.txt"))),
                     P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a1.log2FC_0.1.P30.PBS.TDTvsMIG.TDTup.txt"))))

#
stat.s2.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsCTR.CTRup.txt"))),
                     P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsCTR.TDTup.txt"))),
                     P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsMIG.MIGup.txt"))),
                     P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P06.TDTvsMIG.TDTup.txt"))),
                     P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsCTR.CTRup.txt"))),
                     P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsCTR.TDTup.txt"))),
                     P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsMIG.MIGup.txt"))),
                     P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a2.log2FC_0.2.P30.PBS.TDTvsMIG.TDTup.txt"))))

#
stat.s3.list <- list(P06.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsCTR.CTRup.txt"))),
                     P06.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsCTR.TDTup.txt"))),
                     P06.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsMIG.MIGup.txt"))),
                     P06.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P06.TDTvsMIG.TDTup.txt"))),
                     P30.PBS.TDTvsCTR.CTRup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsCTR.CTRup.txt"))),
                     P30.PBS.TDTvsCTR.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsCTR.TDTup.txt"))),
                     P30.PBS.TDTvsMIG.MIGup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsMIG.MIGup.txt"))),
                     P30.PBS.TDTvsMIG.TDTup = as.vector(unlist(read.table("../subset1/figure1008/subset1.plot7.DEGlist/a3.log2FC_0.3.P30.PBS.TDTvsMIG.TDTup.txt"))))
lapply(stat.s1.list,length)
## $P06.TDTvsCTR.CTRup
## [1] 174
## 
## $P06.TDTvsCTR.TDTup
## [1] 84
## 
## $P06.TDTvsMIG.MIGup
## [1] 70
## 
## $P06.TDTvsMIG.TDTup
## [1] 61
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 28
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 35
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 21
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 5
lapply(stat.s2.list,length)
## $P06.TDTvsCTR.CTRup
## [1] 62
## 
## $P06.TDTvsCTR.TDTup
## [1] 52
## 
## $P06.TDTvsMIG.MIGup
## [1] 11
## 
## $P06.TDTvsMIG.TDTup
## [1] 31
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 8
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 17
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 9
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 4
lapply(stat.s3.list,length)
## $P06.TDTvsCTR.CTRup
## [1] 14
## 
## $P06.TDTvsCTR.TDTup
## [1] 18
## 
## $P06.TDTvsMIG.MIGup
## [1] 3
## 
## $P06.TDTvsMIG.TDTup
## [1] 14
## 
## $P30.PBS.TDTvsCTR.CTRup
## [1] 2
## 
## $P30.PBS.TDTvsCTR.TDTup
## [1] 7
## 
## $P30.PBS.TDTvsMIG.MIGup
## [1] 3
## 
## $P30.PBS.TDTvsMIG.TDTup
## [1] 1

plot8

DEG volcano, set same x/y-axis

pp.DEGs.a$P30.LPS.TDTvsCTR <- DEGs.a.combine$P30.LPS.TDTvsCTR %>%
  mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTR"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.LPS, TDT vs CTR") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")  + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.LPS.TDTvsCTR

pp.DEGs.a$P30.LPS.TDTvsMIG <- DEGs.a.combine$P30.LPS.TDTvsMIG %>%
  mutate(label=ifelse(((p_val_adj < 1e-8 & avg_log2FC<0) | (p_val_adj < 1e-8 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("MIG"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.LPS, TDT vs MIG") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted") + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.LPS.TDTvsMIG

pp.DEGs.a$P30.PBS.TDTvsCTR <- DEGs.a.combine$P30.PBS.TDTvsCTR %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"CTR",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("CTR"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.PBS, TDT vs CTR") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")  + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.PBS.TDTvsCTR

pp.DEGs.a$P30.PBS.TDTvsMIG <- DEGs.a.combine$P30.PBS.TDTvsMIG %>%
  mutate(label=ifelse(((p_val_adj < 1e-3 & avg_log2FC<0) | (p_val_adj < 1e-3 & avg_log2FC>0) | (p_val_adj < 0.05 & abs(avg_log2FC) > 0.5)) &
                        !grepl("^Gm|^Hsp|^Rps|^Rpl|Rik$|Xist|Tsix|Uty|Ddx3y|Eif2s3y|Kdm5d",gene),gene,""),
         sig.=ifelse(p_val_adj<0.05 & avg_log2FC<0,"MIG",ifelse(p_val_adj<0.05 & avg_log2FC>0,"TDT","None")),
         padj=ifelse(p_val_adj<1e-60,p_val_adj+1e-60,p_val_adj)) %>%
  ggplot(aes(y= -log10(padj), x=avg_log2FC, size=pct.1, label=label,fill=sig.)) +
  geom_point(shape=21, alpha=0.65) +
  geom_text_repel(size=2.5, max.overlaps = 500) +
  scale_fill_manual(values = c("MIG"="Skyblue",
                               "TDT"="pink",
                               "None"="grey")) +
  theme_classic() + labs(title="P30.PBS, TDT vs MIG") +
  guides(size = guide_legend(order = 2), fill=guide_legend(order = 1, reverse = F)) +
  geom_vline(xintercept = c(-0.1,0.1), linetype="dotted")  + xlim(c(-1.5,3)) + ylim(c(0,60))
pp.DEGs.a$P30.PBS.TDTvsMIG

plot9

‘LPS_Tdt vs LPS_MIG’, how many DEGs overlap with ‘P06 Tdt vs P06 CTRL’ and ‘P06 Tdt vs MIG’ ?

library(UpSetR)
## Warning: package 'UpSetR' was built under R version 4.0.5
upset(fromList(c(stat.a1.list[3:4],
                 stat.s1.list[1:4])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
  
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))

upset(fromList(c(stat.a1.list[3:4],
                 stat.s1.list[1:2])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
  
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))

upset(fromList(c(stat.a1.list[3:4],
                 stat.s1.list[3:4])),
      order.by = 'freq', nsets = 6, point.size=3.5, line.size =1, text.scale = 2)
  
grid::grid.text("pct.a>0.05, padj<0.05, |log2FC|>0.1",x = 0.75, y=0.95, gp=grid::gpar(fontsize=20))

overlap list

plot10

cut.padj = 0.05
cut.log2FC = 0.3   
  
cut.pct1 = 0.05

df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame()
##   condition cluster up.DEGs
## 1   P30.CTR P30.PBS     425
## 2   P30.CTR P30.LPS     512
## 3   P30.MIG P30.PBS     321
## 4   P30.MIG P30.LPS     391
## 5   P30.TDT P30.PBS     388
## 6   P30.TDT P30.LPS     423
pp.stat.DEG[["b"]] <- df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>%
  group_by(condition,cluster) %>%
  summarise(up.DEGs = n()) %>% as.data.frame() %>%
  ggplot(aes(x=condition, y=up.DEGs, color = cluster)) +
  geom_bar(stat="summary", fun="mean", position = position_dodge(0.75), width = 0.58, fill="white") +
  geom_text(aes(label = up.DEGs),vjust=-0.21, show.legend = F, position = position_dodge(0.75) ) +
  theme_classic(base_size = 15) +
  scale_color_manual(values = c("skyblue","pink"), name="") +
  labs(title=paste0("up.DEGs stat, pct.1>",cut.pct1,", padj<",cut.padj,", |log2FC|>",cut.log2FC), y = "Number of DEGs") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = 9.6),
        title =element_text(size=12, face='bold'))
pp.stat.DEG[["b"]]

sig.LPS <- df.DEGs.b %>% filter(p_val_adj < cut.padj & 
                       abs(avg_log2FC) > cut.log2FC & 
                       pct.1 > cut.pct1) %>% 
  filter(condition == "P30.CTR" & cluster=="P30.LPS")
dim(sig.LPS)
## [1] 512   8
sig.LPS[1:6,]
##          p_val avg_log2FC pct.1 pct.2 p_val_adj cluster     gene condition
## Ccl12        0   3.787020 0.993 0.731         0 P30.LPS    Ccl12   P30.CTR
## Ifitm3       0   2.120287 0.766 0.146         0 P30.LPS   Ifitm3   P30.CTR
## Slfn2        0   1.709130 0.842 0.190         0 P30.LPS    Slfn2   P30.CTR
## Ifi30        0   1.647601 0.966 0.558         0 P30.LPS    Ifi30   P30.CTR
## Hsp90ab1     0   1.617709 1.000 0.937         0 P30.LPS Hsp90ab1   P30.CTR
## Srgn         0   1.541016 0.942 0.499         0 P30.LPS     Srgn   P30.CTR
GEX.seur <- add_geneset_score(GEX.seur, sig.LPS$gene, "LPS.sig_score")
## Summarizing data
pp.feat.10a1 <- FeaturePlot(GEX.seur,features = c("LPS.sig_score"), ncol = 1,
                           raster = T, pt.size = 3.5#,keep.scale = "all"
                           )
pp.feat.10a1

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
pp.feat.10a2 <- FeaturePlot(GEX.seur,features = c("LPS.sig_score"), ncol = 1,
                           raster = T, pt.size = 3.5#,keep.scale = "all"
                           ) & scale_color_gradientn(colors = rev(mapal))
pp.feat.10a2

pp.feat.10b1 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) + geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02) + NoLegend()
pp.feat.10b1

pp.feat.10b2 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"), 
                      same.y.lims = T,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) + NoLegend() & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) &  ylim(c(0,1.1)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.CTR","P30.PBS.CTR")),
                               label.y = c(0.73,0.78,0.87,0.55,0.65,0.75,0.93),
                               size=2.35
                               )
pp.feat.10b2

pp.feat.10c1 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"), 
                      group.by = "new_clusters", #cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) + NoLegend()& geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02) 
pp.feat.10c1

pp.feat.10c2 <- VlnPlot(GEX.seur, features = c("LPS.sig_score"), 
                      group.by = "new_clusters", #cols = color.cnt,
                      ncol = 1, pt.size = 0, raster = F) + NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55)&
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) 
pp.feat.10c2

pp.feat.10c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                          features = c("LPS.sig_score"), 
                        group.by = "cluster_cnt", cols = rep(color.cnt[1:3],3),
                      ncol = 1, pt.size = 0, raster = F)+ NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,0.7)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.a,
                               label.y = (plot.labely.a-1)/10.5,
                               size=1.75
                               ) 
pp.feat.10c3a

pp.feat.10c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                          features = c("LPS.sig_score"), 
                        group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
                      ncol = 1, pt.size = 0, raster = F)  + NoLegend() & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0.35,1)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.b,
                               label.y = (plot.labely.b+1)/8.6,
                               size=1.5
                               )
pp.feat.10c3b

plot11

DAM score

DAM.sig <- read.csv("../../../20230811_10x_LYN/analysis_plus_exogene/figures1002/new/ranking_of_DAM_indicator_genes.csv")
DAM.sig[1:8,]
##   锘縍anking     Gene Frequence Total.score
## 1          1     Spp1        11    35.98221
## 2          2     Apoe        11    31.44704
## 3          3   Ifitm3         9    20.82901
## 4          4 Ifi27l2a         9    20.44047
## 5          5     Lyz2        11    20.22463
## 6          6     Cst7         9    19.48372
## 7          7     Cd74         7    18.24720
## 8          8   Lgals3         7    18.24180
DAM.list <- list(top50=DAM.sig$Gene[1:50],
                 top100=DAM.sig$Gene[1:100],
                 top250=DAM.sig$Gene[1:250],
                 top500=DAM.sig$Gene[1:500])
DAM.list
## $top50
##  [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##  [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
## [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
## [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
## [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
## [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
## [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
## [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
## [49] "Ifit3"    "Apoc1"   
## 
## $top100
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"   
## 
## $top250
##   [1] "Spp1"     "Apoe"     "Ifitm3"   "Ifi27l2a" "Lyz2"     "Cst7"    
##   [7] "Cd74"     "Lgals3"   "Lpl"      "Cd52"     "Ccl5"     "Ccl12"   
##  [13] "Lgals3bp" "Cd63"     "H2-Ab1"   "Fn1"      "H2-Aa"    "Fxyd5"   
##  [19] "Ccl4"     "Cybb"     "Bst2"     "Cstb"     "H2-Eb1"   "Fth1"    
##  [25] "Vim"      "Tspo"     "Ctsb"     "Ccl3"     "Axl"      "Anxa5"   
##  [31] "Isg15"    "Lgals1"   "Ccl2"     "Ifi204"   "Igf1"     "Ch25h"   
##  [37] "Mif"      "Cxcl10"   "Plin2"    "Fabp5"    "Il1b"     "Crip1"   
##  [43] "Slfn2"    "B2m"      "Irf7"     "Cd72"     "Capg"     "Ms4a6c"  
##  [49] "Ifit3"    "Apoc1"    "Fcgr4"    "Il1rn"    "Wfdc17"   "Cxcl2"   
##  [55] "Cxcl16"   "Prdx1"    "Rtp4"     "H2-D1"    "Pkm"      "Stat1"   
##  [61] "Anxa2"    "Gpnmb"    "Zbp1"     "Ftl1"     "Ldha"     "Npc2"    
##  [67] "Ly6a"     "Oasl2"    "Gapdh"    "Prdx5"    "Gbp2"     "Grn"     
##  [73] "Ifi207"   "Ifitm2"   "Tlr2"     "Txn1"     "Phf11b"   "Ctsz"    
##  [79] "Pianp"    "Cd36"     "Itgax"    "Fgl2"     "Ccl6"     "Iigp1"   
##  [85] "Ccl7"     "H2-K1"    "Pld3"     "Tpi1"     "Akr1a1"   "Usp18"   
##  [91] "Rab7b"    "Tmsb10"   "Cd44"     "Aldoa"    "Hcar2"    "Acod1"   
##  [97] "Cd300lb"  "Ctsc"     "Gpr65"    "H2-Q7"    "Cdkn1a"   "Psat1"   
## [103] "Trim30a"  "Cxcl14"   "Srgn"     "Mmp12"    "Bcl2a1b"  "Tap1"    
## [109] "AB124611" "Xaf1"     "Ly6e"     "Psme1"    "Ctsl"     "Sp100"   
## [115] "Cxcr4"    "Psmb8"    "AA467197" "Pgk1"     "Emp3"     "Csf1"    
## [121] "Mcemp1"   "Gusb"     "Pim1"     "Ctse"     "Cox4i1"   "Il12b"   
## [127] "Msrb1"    "Npm1"     "Alcam"    "Psme2"    "Apoc2"    "Bhlhe40" 
## [133] "Stmn1"    "Gm4951"   "Pla2g7"   "Plaur"    "Tor3a"    "Hspe1"   
## [139] "Tpt1"     "Ndufa1"   "Flt1"     "Tgfbi"    "Cox6c"    "Irgm1"   
## [145] "Ifit2"    "Uba52"    "Igf2r"    "Bola2"    "Ank"      "Tyrobp"  
## [151] "Tpm4"     "Ass1"     "Ms4a4c"   "Ifit1"    "Ybx1"     "Sod2"    
## [157] "Cox8a"    "Fam20c"   "Oas1a"    "Arg1"     "Ms4a7"    "Smpdl3a" 
## [163] "Sh3pxd2b" "Fau"      "Gnas"     "Phf11d"   "Ehd1"     "Saa3"    
## [169] "Cox5a"    "Atox1"    "Id3"      "Lrpap1"   "Olr1"     "Sh3bgrl3"
## [175] "Sash1"    "Hint1"    "Il2rg"    "Ctsd"     "Postn"    "H2-T23"  
## [181] "Ucp2"     "Sdc3"     "Uqcrq"    "Cox6a1"   "Cd300lf"  "Syngr1"  
## [187] "Timp2"    "Atp5e"    "Id2"      "S100a6"   "Cd14"     "Tubb6"   
## [193] "Anp32b"   "Fcgr2b"   "Cd83"     "Psmb9"    "Bcl2a1a"  "Aprt"    
## [199] "Mfsd12"   "Psap"     "Cox6b1"   "Lilr4b"   "Plac8"    "Glrx"    
## [205] "Scimp"    "Rilpl2"   "Psmb6"    "Gm11808"  "Chmp4b"   "Actr3b"  
## [211] "Ly86"     "Fundc2"   "Ifi211"   "Hif1a"    "Serf2"    "Dram1"   
## [217] "Parp14"   "Ptgs2"    "Cxcl9"    "Myo5a"    "Gabarap"  "Cyp4f18" 
## [223] "Shisa5"   "Lilrb4a"  "Cpd"      "Iqgap1"   "Slfn5"    "Tnfaip2" 
## [229] "Nme1"     "Cotl1"    "Tagln2"   "Mt1"      "Mpeg1"    "C3"      
## [235] "Ube2l6"   "Clic4"    "Naaa"     "Gas7"     "Sdcbp"    "Nampt"   
## [241] "Ell2"     "Samhd1"   "Rtcb"     "Eef1g"    "Hmgb2"    "Gng5"    
## [247] "Nfil3"    "Adora1"   "Tpd52"    "Ptger4"  
## 
## $top500
##   [1] "Spp1"          "Apoe"          "Ifitm3"        "Ifi27l2a"     
##   [5] "Lyz2"          "Cst7"          "Cd74"          "Lgals3"       
##   [9] "Lpl"           "Cd52"          "Ccl5"          "Ccl12"        
##  [13] "Lgals3bp"      "Cd63"          "H2-Ab1"        "Fn1"          
##  [17] "H2-Aa"         "Fxyd5"         "Ccl4"          "Cybb"         
##  [21] "Bst2"          "Cstb"          "H2-Eb1"        "Fth1"         
##  [25] "Vim"           "Tspo"          "Ctsb"          "Ccl3"         
##  [29] "Axl"           "Anxa5"         "Isg15"         "Lgals1"       
##  [33] "Ccl2"          "Ifi204"        "Igf1"          "Ch25h"        
##  [37] "Mif"           "Cxcl10"        "Plin2"         "Fabp5"        
##  [41] "Il1b"          "Crip1"         "Slfn2"         "B2m"          
##  [45] "Irf7"          "Cd72"          "Capg"          "Ms4a6c"       
##  [49] "Ifit3"         "Apoc1"         "Fcgr4"         "Il1rn"        
##  [53] "Wfdc17"        "Cxcl2"         "Cxcl16"        "Prdx1"        
##  [57] "Rtp4"          "H2-D1"         "Pkm"           "Stat1"        
##  [61] "Anxa2"         "Gpnmb"         "Zbp1"          "Ftl1"         
##  [65] "Ldha"          "Npc2"          "Ly6a"          "Oasl2"        
##  [69] "Gapdh"         "Prdx5"         "Gbp2"          "Grn"          
##  [73] "Ifi207"        "Ifitm2"        "Tlr2"          "Txn1"         
##  [77] "Phf11b"        "Ctsz"          "Pianp"         "Cd36"         
##  [81] "Itgax"         "Fgl2"          "Ccl6"          "Iigp1"        
##  [85] "Ccl7"          "H2-K1"         "Pld3"          "Tpi1"         
##  [89] "Akr1a1"        "Usp18"         "Rab7b"         "Tmsb10"       
##  [93] "Cd44"          "Aldoa"         "Hcar2"         "Acod1"        
##  [97] "Cd300lb"       "Ctsc"          "Gpr65"         "H2-Q7"        
## [101] "Cdkn1a"        "Psat1"         "Trim30a"       "Cxcl14"       
## [105] "Srgn"          "Mmp12"         "Bcl2a1b"       "Tap1"         
## [109] "AB124611"      "Xaf1"          "Ly6e"          "Psme1"        
## [113] "Ctsl"          "Sp100"         "Cxcr4"         "Psmb8"        
## [117] "AA467197"      "Pgk1"          "Emp3"          "Csf1"         
## [121] "Mcemp1"        "Gusb"          "Pim1"          "Ctse"         
## [125] "Cox4i1"        "Il12b"         "Msrb1"         "Npm1"         
## [129] "Alcam"         "Psme2"         "Apoc2"         "Bhlhe40"      
## [133] "Stmn1"         "Gm4951"        "Pla2g7"        "Plaur"        
## [137] "Tor3a"         "Hspe1"         "Tpt1"          "Ndufa1"       
## [141] "Flt1"          "Tgfbi"         "Cox6c"         "Irgm1"        
## [145] "Ifit2"         "Uba52"         "Igf2r"         "Bola2"        
## [149] "Ank"           "Tyrobp"        "Tpm4"          "Ass1"         
## [153] "Ms4a4c"        "Ifit1"         "Ybx1"          "Sod2"         
## [157] "Cox8a"         "Fam20c"        "Oas1a"         "Arg1"         
## [161] "Ms4a7"         "Smpdl3a"       "Sh3pxd2b"      "Fau"          
## [165] "Gnas"          "Phf11d"        "Ehd1"          "Saa3"         
## [169] "Cox5a"         "Atox1"         "Id3"           "Lrpap1"       
## [173] "Olr1"          "Sh3bgrl3"      "Sash1"         "Hint1"        
## [177] "Il2rg"         "Ctsd"          "Postn"         "H2-T23"       
## [181] "Ucp2"          "Sdc3"          "Uqcrq"         "Cox6a1"       
## [185] "Cd300lf"       "Syngr1"        "Timp2"         "Atp5e"        
## [189] "Id2"           "S100a6"        "Cd14"          "Tubb6"        
## [193] "Anp32b"        "Fcgr2b"        "Cd83"          "Psmb9"        
## [197] "Bcl2a1a"       "Aprt"          "Mfsd12"        "Psap"         
## [201] "Cox6b1"        "Lilr4b"        "Plac8"         "Glrx"         
## [205] "Scimp"         "Rilpl2"        "Psmb6"         "Gm11808"      
## [209] "Chmp4b"        "Actr3b"        "Ly86"          "Fundc2"       
## [213] "Ifi211"        "Hif1a"         "Serf2"         "Dram1"        
## [217] "Parp14"        "Ptgs2"         "Cxcl9"         "Myo5a"        
## [221] "Gabarap"       "Cyp4f18"       "Shisa5"        "Lilrb4a"      
## [225] "Cpd"           "Iqgap1"        "Slfn5"         "Tnfaip2"      
## [229] "Nme1"          "Cotl1"         "Tagln2"        "Mt1"          
## [233] "Mpeg1"         "C3"            "Ube2l6"        "Clic4"        
## [237] "Naaa"          "Gas7"          "Sdcbp"         "Nampt"        
## [241] "Ell2"          "Samhd1"        "Rtcb"          "Eef1g"        
## [245] "Hmgb2"         "Gng5"          "Nfil3"         "Adora1"       
## [249] "Tpd52"         "Ptger4"        "Eif2ak2"       "Areg"         
## [253] "Rsad2"         "Slc31a1"       "Gm2000"        "Cox7c"        
## [257] "Tmem256"       "Cox5b"         "Cyba"          "Il18bp"       
## [261] "Selenow"       "Myl6"          "H2-DMa"        "Rai14"        
## [265] "Dab2"          "Hmox1"         "Tmsb4x"        "Ndufa2"       
## [269] "Cd68"          "Ccdc86"        "Atp6v1a"       "Cox7a2"       
## [273] "Calm1"         "Uqcrh"         "Socs2"         "Gpi1"         
## [277] "Ubl5"          "Colec12"       "Ifit3b"        "Scpep1"       
## [281] "S100a11"       "F3"            "Ms4a6d"        "Hacd2"        
## [285] "Chil3"         "Tuba1c"        "Rap2b"         "Gp49a"        
## [289] "Milr1"         "Fcgr1"         "Stat2"         "Atp5j2"       
## [293] "Uqcr10"        "Ssr4"          "Bcar3"         "Gm49368"      
## [297] "Tmem106a"      "Tubb5"         "Naca"          "Hspa8"        
## [301] "Atp5k"         "Bag1"          "Sec61b"        "Gyg"          
## [305] "Cox7b"         "Ly6c2"         "Top2a"         "Aldh2"        
## [309] "Dpp7"          "Eif3k"         "Cd48"          "C4b"          
## [313] "Pycard"        "Atp5l"         "Uqcr11"        "Osm"          
## [317] "Prelid1"       "Rnf213"        "Prdx4"         "Arpc1b"       
## [321] "Ndufv3"        "Sp110"         "Ndufa13"       "Abca1"        
## [325] "Gm1673"        "Wdr89"         "Sh3bgrl"       "Smdt1"        
## [329] "Gatm"          "Gpr84"         "Slamf8"        "Ccrl2"        
## [333] "Tomm7"         "Gas8"          "Ly6i"          "Bcl2a1d"      
## [337] "Esd"           "Dhx58"         "Ctsa"          "Rxrg"         
## [341] "Prdx6"         "Ndufc1"        "Polr2f"        "Sdc4"         
## [345] "Atp5j"         "Litaf"         "Atp6v0e"       "Cspg4"        
## [349] "Ranbp1"        "Ifi35"         "Fcer1g"        "Calm3"        
## [353] "Rhoc"          "Pde4b"         "Atp5g1"        "Gpx3"         
## [357] "Psmb1"         "Gpx1"          "Eef1a1"        "Ly9"          
## [361] "Igtp"          "H2-Q6"         "Herc6"         "Cd9"          
## [365] "Ran"           "Hebp1"         "Eno1"          "Cdk18"        
## [369] "Eef1b2"        "Gbp7"          "Glipr1"        "Atp6v1f"      
## [373] "H2-DMb1"       "Btf3"          "Slc25a3"       "Brdt"         
## [377] "Csf2ra"        "Eif3f"         "Hpse"          "Gm14305"      
## [381] "Gm14410"       "H2-Q4"         "Ndufb9"        "Epsti1"       
## [385] "Dnaaf3"        "Pf4"           "S100a4"        "Atp5h"        
## [389] "Apobec3"       "Hsp90ab1"      "Tnf"           "Ctss"         
## [393] "Gas6"          "Gbp3"          "Gng12"         "Nceh1"        
## [397] "Mgst1"         "Cfl1"          "Dtnbp1"        "Slc2a6"       
## [401] "Eif5a"         "Atp5c1"        "ptchd1"        "Ptchd1"       
## [405] "Psma7"         "Lap3"          "Rbm3"          "Cycs"         
## [409] "Cox6a2"        "Abcg1"         "Prr15"         "Ahnak"        
## [413] "Ndufb7"        "Nrp1"          "Lmnb1"         "Ninj1"        
## [417] "Mrpl33"        "Arpc3"         "Phlda1"        "Acp5"         
## [421] "Atp5g3"        "C5ar1"         "Arl5c"         "Parp9"        
## [425] "Ndufb8"        "Txndc17"       "Tbca"          "Gm49339"      
## [429] "Tma7"          "Fblim1"        "Eif3h"         "Psme2b"       
## [433] "Mrpl30"        "Arpc2"         "Ptma"          "Rac2"         
## [437] "Ctsh"          "Dcstamp"       "Clec4n"        "Dbi"          
## [441] "H2-T22"        "Sem1"          "Msr1"          "Bax"          
## [445] "S100a10"       "Fabp3"         "Snrpg"         "Bcl2a1"       
## [449] "Serpine2"      "Snrpd2"        "Cd180"         "Pgam1"        
## [453] "2310061I04Rik" "S100a1"        "Serpinb6a"     "Actg1"        
## [457] "Uqcc2"         "Tuba1b"        "Neurl2"        "Gcnt2"        
## [461] "Smc2"          "Atp5b"         "Ifih1"         "Nap1l1"       
## [465] "Cd274"         "Rgs1"          "Parp12"        "Serpine1"     
## [469] "Nsa2"          "Cebpb"         "Csf2rb"        "Cfb"          
## [473] "Ndufa6"        "Sdf2l1"        "Anxa4"         "Psma5"        
## [477] "Nfkbiz"        "Ctla2a"        "Atp5o"         "Ube2l3"       
## [481] "Lipa"          "Pfn1"          "Ndufa7"        "Ndufs6"       
## [485] "Psmb10"        "Mapkapk2"      "Aif1"          "Smc4"         
## [489] "Itgb1bp1"      "Nme2"          "Pfdn5"         "Rassf3"       
## [493] "1810058I24Rik" "Pgls"          "Clec4d"        "Mrpl54"       
## [497] "Tmem160"       "Pomp"          "Slc7a11"       "F10"
lapply(DAM.list, length)
## $top50
## [1] 50
## 
## $top100
## [1] 100
## 
## $top250
## [1] 250
## 
## $top500
## [1] 500
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top50, "DAM.sig_top50")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top100, "DAM.sig_top100")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top250, "DAM.sig_top250")
## Summarizing data
GEX.seur <- add_geneset_score(GEX.seur, DAM.list$top500, "DAM.sig_top500")
## Summarizing data
pp.feat.11a1 <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                           raster = T, pt.size = 3.5#,keep.scale = "all"
                           )
pp.feat.11a1

mapal <- colorRampPalette(RColorBrewer::brewer.pal(4,"Spectral"))(120)
pp.feat.11a2 <- FeaturePlot(GEX.seur,features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), ncol = 4,
                        raster = T, pt.size = 3.5,
            keep.scale = "all") & scale_color_gradientn(colors = rev(mapal))
pp.feat.11a2

pp.feat.11b1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.feat.11b1

pp.feat.11b2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      same.y.lims = T,
                      group.by = "cnt", cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & 
    geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55) & 
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.TDT","P30.PBS.TDT")),
                               label.y = c(0.6,0.75,0.9,0.6,0.75,0.9,1.05),
                               size=2.35
                               )
pp.feat.11b2

pp.feat.11c1 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      group.by = "new_clusters", #cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_jitter(alpha=0.15, shape=16, width = 0.2, size = 0.02)
pp.feat.11c1

pp.feat.11c2 <- VlnPlot(GEX.seur, features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                      group.by = "new_clusters", #cols = color.cnt,
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.2, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=3, color="black", alpha=0.55)
pp.feat.11c2

pp.feat.11c3a <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.PBS") & new_clusters %in% paste0("C",1:3)), 
                          features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                        group.by = "cluster_cnt", cols = rep(color.cnt[1:3],3),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(-0.4,0.65)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.a,
                               label.y = (plot.labely.a-1)/10.5,
                               size=1.75
                               )
pp.feat.11c3a

pp.feat.11c3b <- VlnPlot(subset(GEX.seur,subset= age %in% c("P30.LPS") & new_clusters %in% paste0("C",4:11)), 
                          features = c("DAM.sig_top50","DAM.sig_top100","DAM.sig_top250","DAM.sig_top500"), 
                        group.by = "cluster_cnt", cols = rep(color.cnt[4:6],8),
                      ncol = 2, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(-0.25,1.3)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = plot.contrast.b,
                               label.y = (plot.labely.b-1)/5.5,
                               size=1.5
                               ) 
pp.feat.11c3b

#saveRDS(GEX.seur,"./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")
#GEX.seur <- readRDS("./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")

newlist1009

dotplot and violin in conditions

GEX.seur
## An object of class Seurat 
## 17900 features across 15536 samples within 1 assay 
## Active assay: RNA (17900 features, 1437 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
marker1009 <- as.vector(unlist(read.table("./dotplot.subset2_newlist1009.txt")))
marker1009 <- unique(marker1009)
marker1009
##  [1] "Ms4a6c"  "Fn1"     "Apoe"    "Ms4a6d"  "Spp1"    "Ms4a4c"  "Lilrb4a"
##  [8] "C4b"     "Tspo"    "Ms4a6b"  "Axl"     "Lpl"     "Cxcl16"  "Ccl2"   
## [15] "Ccl12"   "Saa3"    "Ccl5"    "Il1b"    "Ccl3"    "Cxcl10"  "Ccl4"   
## [22] "Il1a"    "Tnf"     "Fcrls"   "Gpr34"   "Siglech" "P2ry12"  "P2ry13" 
## [29] "Tmem119" "Sall1"   "Slc2a5"  "Cst3"
length(marker1009)
## [1] 32
sum(marker1009 %in% markers.new$gene)
## [1] 31
marker1009.sort <- (markers.new %>% group_by(cluster) %>% 
                  filter(gene %in% marker1009) %>%
                   ungroup() %>%
    #arrange(desc(avg_log2FC*pct.1),gene) %>%
    arrange(desc(avg_log2FC),gene) %>%
                             distinct(gene, .keep_all = TRUE) %>%
                             arrange(cluster,p_val_adj))$gene
length(marker1009)
## [1] 32
length(marker1009.sort)
## [1] 31
setdiff(marker1009,marker1009.sort)
## [1] "Fcrls"
marker1009.new <- c(marker1009.sort[1:8],"Fcrls",marker1009.sort[9:31])
pp.dot.x1a <- DotPlot(GEX.seur, features = marker1009, group.by = "new_clusters"
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.x1a

pp.dot.x1a <- DotPlot(GEX.seur, features = marker1009.new, group.by = "new_clusters",cols = c("lightgrey","red")
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Clusters")
pp.dot.x1a

pp.dot.x1a  + scale_color_gradientn(colours  = rev(mapal))

pp.dot.x1b <- DotPlot(GEX.seur, features = marker1009.new, group.by = "cnt",cols = c("lightgrey","red")
                         ) +
                 theme(axis.text.x = element_text(angle = 60, hjust = 1, vjust = 1)) + labs(y="Conditions")
pp.dot.x1b

pp.dot.x1b  + scale_color_gradientn(colours  = rev(mapal))

pp.violin.x1b <- VlnPlot(GEX.seur, features = marker1009.new, group.by = "cnt", cols = color.cnt,
                      ncol = 5, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55)
pp.violin.x1b

pp.violin.x1c <- VlnPlot(GEX.seur, features = marker1009.new, group.by = "cnt", cols = color.cnt,
                      ncol = 5, pt.size = 0, raster = F) & geom_boxplot(outlier.size = 0, fill="white", width=0.3, size=0.1, alpha=0.55) &
  stat_summary(fun=mean, geom="point", shape=18, size=2, color="black", alpha=0.55) & ylim(c(0,7.2)) &
  ggpubr::stat_compare_means(aes(lable = ..p.signif..), 
                               method = "wilcox.test",
                               comparisons = list(c("P30.LPS.CTR","P30.LPS.MIG"),
                                                  c("P30.LPS.MIG","P30.LPS.TDT"),
                                                  c("P30.LPS.CTR","P30.LPS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.MIG"),
                                                  c("P30.PBS.MIG","P30.PBS.TDT"),
                                                  c("P30.PBS.CTR","P30.PBS.TDT"),
                                                  c("P30.LPS.TDT","P30.PBS.TDT")),
                               label.y = c(6.3,6.8,7.3,6.3,6.8,7.3,7.8)-1,
                               size=1.75
                               )
pp.violin.x1c

forGEO

GEX.seur <- readRDS("./Spp1tdt.subset2_P30LPS_P30PBS..new1008.rds")
GEX.seur
## An object of class Seurat 
## 17900 features across 15536 samples within 1 assay 
## Active assay: RNA (17900 features, 1437 variable features)
##  3 dimensional reductions calculated: pca, tsne, umap
GEX.seur@meta.data[,grep("snn|pANN",colnames(GEX.seur@meta.data),value = T)] <- NULL
head(GEX.seur@meta.data)
##                    orig.ident nCount_RNA nFeature_RNA percent.mt percent.rb
## AAACCCAAGGTCGACA-1    Spp1tdt       5001         2062   2.699460  12.437512
## AAACCCAAGGTTCTTG-1    Spp1tdt       6501         2457   2.553453  12.305799
## AAACCCACAAACTCGT-1    Spp1tdt       4188         1841   2.960840   8.572111
## AAACCCACACGCGGTT-1    Spp1tdt       4524         1913   1.215738   9.858532
## AAACCCACATCCGTGG-1    Spp1tdt       2337         1225   3.979461   7.659392
## AAACCCAGTATCGAAA-1    Spp1tdt       4318         1955   2.292728  11.741547
##                         FB.info       S.Score   G2M.Score Phase seurat_clusters
## AAACCCAAGGTCGACA-1  P30.LPS.CTR  0.0061184939 -0.06523366     S               0
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT1  0.0002076412 -0.07102627     S               0
## AAACCCACAAACTCGT-1  P30.PBS.MIG  0.0192137320  0.06031440   G2M               7
## AAACCCACACGCGGTT-1 P30.LPS.TDT2 -0.0535299003 -0.05464656    G1               0
## AAACCCACATCCGTGG-1  P30.PBS.TDT  0.0489756368 -0.04211671     S               2
## AAACCCAGTATCGAAA-1  P30.LPS.MIG -0.0314645626  0.02183367   G2M               0
##                    DoubletFinder0.05 DoubletFinder0.1 sort_clusters     age
## AAACCCAAGGTCGACA-1           Singlet          Singlet             0 P30.LPS
## AAACCCAAGGTTCTTG-1           Singlet          Singlet             0 P30.LPS
## AAACCCACAAACTCGT-1           Singlet          Singlet             7 P30.PBS
## AAACCCACACGCGGTT-1           Singlet          Singlet             0 P30.LPS
## AAACCCACATCCGTGG-1           Singlet          Singlet             2 P30.PBS
## AAACCCAGTATCGAAA-1           Singlet          Singlet             0 P30.LPS
##                            cnt   preAnno       FB.new LPS.TDT.up  LPS.TDT.dn
## AAACCCAAGGTCGACA-1 P30.LPS.CTR Microglia  P30.LPS.CTR  0.2671683  0.24670044
## AAACCCAAGGTTCTTG-1 P30.LPS.TDT Microglia P30.LPS.TDT1  0.2190005  0.21501295
## AAACCCACAAACTCGT-1 P30.PBS.MIG Microglia  P30.PBS.MIG  0.3136800  0.02545367
## AAACCCACACGCGGTT-1 P30.LPS.TDT Microglia P30.LPS.TDT2  0.2974201  0.09338095
## AAACCCACATCCGTGG-1 P30.PBS.TDT Microglia  P30.PBS.TDT  0.2758054 -0.03406151
## AAACCCAGTATCGAAA-1 P30.LPS.MIG Microglia  P30.LPS.MIG  0.2300058  0.11176488
##                    new_clusters    age1    cluster_cnt DAM.sig_top50
## AAACCCAAGGTCGACA-1           C7 P30.LPS C7_P30.LPS.CTR    0.12647123
## AAACCCAAGGTTCTTG-1           C7 P30.LPS C7_P30.LPS.TDT    0.04136910
## AAACCCACAAACTCGT-1           C2 P30.PBS C2_P30.PBS.MIG   -0.03826143
## AAACCCACACGCGGTT-1           C7 P30.LPS C7_P30.LPS.TDT    0.17479143
## AAACCCACATCCGTGG-1           C2 P30.PBS C2_P30.PBS.TDT    0.06229684
## AAACCCAGTATCGAAA-1           C7 P30.LPS C7_P30.LPS.MIG    0.19221428
##                    DAM.sig_top100 DAM.sig_top250 DAM.sig_top500 LPS.sig_score
## AAACCCAAGGTCGACA-1     0.09610184     0.09546331      0.2497200     0.5988063
## AAACCCAAGGTTCTTG-1     0.02314892     0.06780191      0.2634468     0.6298800
## AAACCCACAAACTCGT-1    -0.02304032    -0.01821938      0.1498645     0.2308747
## AAACCCACACGCGGTT-1     0.26609309     0.19238189      0.3066293     0.6339013
## AAACCCACATCCGTGG-1    -0.06607385    -0.05013960      0.1023302     0.1488521
## AAACCCAGTATCGAAA-1     0.07223828     0.07965503      0.2186862     0.4894010
#saveRDS(GEX.seur,"I:/Shared_win/projects/202310_Spp1DAM/forGEO/seur_obj/Spp1TDT_P30LPS.final.seur_obj.rds")